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ABSTRACT 


The  presence  of  air  bubbles  in  ship  wakes  and  their  dynamic  interaction  with  the 
turbulent  vortical  flows  create  security  problems  by  modifying  the  hydro-acoustic 
properties  of  ship  navigation.  This  study  concerns  understanding  and  controlling  the 
mechanisms,  which  may  lead  to  selective  concentration  of  bubbles  to  form  clusters  or 
clouds  and  to  predict  their  size  distribution  and  motion.  Numerical  simulations  were 
conducted  using  the  Large  Eddy  Simulation  (LES)  technique  in  conjunction  with  a 
Lagrangian  particle  tracking  (LPT)  technique  appropriate  for  dispersed  two-phase 
turbulent  flows. 

In  order  to  cut  down  on  the  execution  time  of  LES,  the  simulation  for  the  flow 
around  the  ship  model  are  not  considered,  instead  the  simulations  started  just  after  the 
ship  where  the  inlet  conditions  are  prescribed  with  the  help  of  a  newly  developed 
Random  Flow  Generation  (RFG)  procedure,  and  to  compute  bubble  distributions. 
Moreover,  to  improve  turn-around  time  of  the  computations  and  allow  LES  of  the 
developing  wake  in  a  larger  domain,  parallel  simulation  tools  are  developed  and  adopted 
in  this  study. 

The  turbulence  characteristics  of  ship  wakes  on  a  straight  tract  and  a  circular  tract 
are  investigated  using  the  above  mentioned  techniques,  and  processes  (e.g.  free  surface 
effects,  anisotropy,  etc.)  contributing  to  turbulence  generation  are  identified,  and 
appropriate  sub-grid  scale  (SGS)  models  are  developed  and  applied.  For  the  first  time,  a 
relatively  long  developing  near  wake  of  three  ship  lengths  was  simulated  using  LES  on 
parallel  computers  with  more  than  six  million  nodes. 
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1  INTRODUCTION 


1.1  Background 

The  problem  of  bubble  dynamics  in  ship  wakes  is  important  in  naval  hydrodynamics 
because  of  ship  security  as  some  bubbles  somehow  form  clusters  within  the  turbulent  eddies 
(vortical  structures)  and  some  are  able  to  persist  for  miles  leaving  a  field  signature.  Bubbles  can 
be  generated  by  various  mechanisms,  e.g.  cavitation,  wave  breaking,  ship  nose  hydrodynamics, 
and  propeller  ventilation.  The  acoustic  properties  such  as  the  resonance  frequency  of  the  media 
are  a  function  of  local  bubble  size  and  population  as  well  as  the  random  motion  and  pressure 
fluctuations  of  the  surrounding  liquid.  The  detailed  microscale  dynamics  of  turbulent  bubbly 
flows  leading  to  acoustic  radiation  needs  to  be  studied  in  order  to  understand  and  eventually  to 
control  hydro-acoustics  properties.  Unless  the  vortical  structures  at  the  lower  end  of  the  integral 
length  scale  of  turbulence  are  resolved  in  relation  to  bubble  size  and  cluster  size  the  bubble 
dynamics  can  not  be  understood.  Most  researchers  use  classical  turbulence  models  called 
Reynolds  Averaged  Navier-Stokes  (RANS)  models.  The  classical  turbulence  models  need 
extensive  empirical  input  and  their  success  is  quite  limited.  The  direct  numerical  simulations  can 
provide  very  detailed  information,  but  application  to  practical  problems  such  as  that  of  the  wake 
of  a  ship  is  impossible  due  to  excessive  amount  of  computer  time  and  storage  requirements. 
Large-eddy  simulation  (LES)  on  the  other  hand  resolves  only  the  large  eddies  to  the  extent  it  is 
necessary,  hence  optimizing  details,  accuracy,  and  computational  requirements.  Not  many 
researchers  utilize  large  eddy  simulations  for  this  problem  in  spite  of  the  successful  predictions 
seen  in  other  areas  such  as  atmospheric  turbulence,  flow  past  bluff  bodies,  boundary  and  shear 
layers  (see  Galperin  and  Orszag,  1993;  Lesieur  and  Metais,  1996;  Reynolds,  1989).  The 
interaction  of  two-phases  from  the  perspective  of  turbulence  enhancement  or  damping  in  the  near 
wake  of  ships  seems  to  play  an  active  role  in  determining  bubble  size  and  population.  This,  in 
turn,  determines  the  phenomenon  in  the  far  wake  which  can  be  predicted  using  non-conventional 
techniques  such  as  large  eddy  simulations.  This  is  the  focus  area  of  the  present  study. 
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1.2  Current  Status  of  Turbulent  Two-phase  Flow  Modeling 

We  shall  consider  only  the  non-reacting,  isothermal  two-phase  flows  with  a  dispersed 
phase  within  a  continuous  phase.  The  local  instantaneous  formulation  for  such  flows  is  fairly 
well  established  (Anderson  and  Jackson,  1967;  Ishii,  1975;  Drew,  1992;  Joseph  and  Lundgren, 
1990)  The  continuum  equations  are  identical  to  the  Navier-Stokes  Equations  for  a  variable 
property  fluid  except  for  the  phase  interaction  terms.  The  constitutive  equations  for  the 
interaction  terms  are  not  as  well  understood,  and  there  is  still  much  research  and  debate 
underway  for  the  advancement  of  these  equations.  The  equation  of  motion  for  a  single  particle  is 
the  well  known  BBO  (Basset-Boussinesq-Osceen)  equation  which  has  been  extended  to  the  case 
of  unsteady,  spatially  non-uniform  ambient  flow  (Maxey  and  Riley,  1983;  Maxey,  1987). 
Simplified  versions  and  many  variants  of  this  equation  are  used  in  applications  to  specific 
problems.  For  example,  when  the  fluid  to  particle  density  ratio  is  small,  many  terms  in  the 
particle  equation  can  be  neglected  leading  to  a  fairly  simple  equation  involving  only  the  drag  and 
lift  forces. 

Numerical  solution  of  the  turbulent  single  phase  flow  equations  can  be  categorized  into 
three  major  groups,  (1)  Reynolds  averaged  Navier-Stokes  (RANS)  approach,  (2)  large-eddy 
simulation  (LES)  technique,  and  (3)  direct  numerical  simulation  (DNS)  approach.  The  same 
approaches  have  been  extended  to  two-phase  flows  but  with  much  more  controversy  and 
empiricism  primarily  due  to  the  complex  nature  of  the  interaction  terms.  Before  we  briefly 
explain  the  merits  of  these  three  approaches  another  classification  of  two-phase  flow  simulations 
needs  to  be  mentioned,  these  are:  the  Eulerian-Lagrangian  (EL),  Eulerian-Eulerian,  (EE)  and 
Mixed  Eulerian-Lagrangian  (MEL)  formulations. 

In  the  EL  approach,  the  “Monte-Carlo  Method”  is  used  almost  exclusively  (see  e.g.  Chen 
and  Crowe,  1984;  Mostafaand  Mongia,  1987,  1989;  Ounis  and  Ahmadi,  1990;  Gouesbet,  1992). 
The  equation  of  motion  for  the  particle  is  solved  directly  using  a  randomly  generated  fluid 
velocity  fluctuation  added  to  the  calculated  mean  velocity  field.  The  fluctuating  components  of 
the  fluid  velocity  are  calculated  using  a  probability  density  distribution  function  (PDF)  which  is 
usually  assumed  to  be  Gaussian.  Its  mean  is  equal  to  the  time  averaged  fluid  velocity,  and  the 
variance  is  the  root  mean  square  of  the  fluid  velocity  functions,  rms  ug  =  (2k/3)  ,  where  k  is  the 
turbulent  kinetic  energy.  This  relation  is  only  valid  for  locally  isotropic  turbulence.  The  time 
variation  of  the  fluid  velocity  fluctuations  are  calculated  by  random  sampling  of  the  PDF  at 
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certain  time  intervals.  This  time  interval  is  called  the  particle  fluid  interaction  time,  tjnt.  It  is  taken 
as  the  minimum  of  a  characteristic  eddy  life  time,  te,  and  a  particle  transit  time,  tr,  defined  as  the 
time  it  takes  a  particle  to  traverse  across  an  eddy  of  size,  le,  which  is  usually  related  to  e,  the 
dissipation  rate  of  k. 

In  the  EE  approach  the  motion  of  the  particulate  phase  is  resolved  by  solving  continuum 
transport  equations  (continuity  and  momentum)  similar  to  those  for  the  fluid  phase  in  an  Eulerian 
frame  of  reference  (Durst  et  al.,  1984;  Mostafa  and  Mongia,  1987;  Chen  and  Pereira,  1995).  Here 
the  major  problem  is  the  modeling  of  turbulence  augmentation  by  the  presence  of  the  dispersed 
phase  in  addition  to  the  modeling  of  interfacial  momentum  exchange  between  the  two  phases. 

A  review  and  comparison  of  the  Lagrangian  and  Eulerian  Approaches  is  presented  by 
Durst  et  al.  (1984),  Mostafa  and  Mongia  (1987),  and  Gouesbet  (1992).  Durst  et  al.  concludes  that 
for  the  Lagrangian  approach  the  computer  storage  requirement  does  not  increase  with  the 
number  of  particle  size  groups.  In  order  to  obtain  more  complete  information  in  an  Eulerian 
approach,  the  equation  of  motion  for  each  representative  particle  size  group  must  be  solved.  This 
would  be  very  time  consuming.  Mostafa  and  Mongia  concludes  that  the  Monte-Carlo  technique 
is  more  expensive  than  the  multi-size  Eulerian  treatment.  Elgobashi  (1996)  points  out  that  in 
Lagrangian  handling  of  the  particle  equation  of  motion  about  90%  time  is  spent  in  particle 
tracking.  However,  this  situation  can  be  improved  considerably  if  the  Lagrangian  particle 
equations  are  solved  in  parallel  using  efficient  parallelization  techniques.  A  further  advantage  of 
the  Lagrangian  approach  is  that  droplet  or  bubble  coalescence  and/or  breakup  (disintegration) 
can  be  modeled  with  relative  ease. 

In  the  MEL  approach,  the  equations  of  mass,  momentum,  and  energy  along  with  the 
trajectory  equations  for  a  range  of  particle  sizes  or  injection  locations  and  properties,  are  solved 
in  Lagrangian  frame  of  reference  to  determine  the  history  of  each  particle.  The  mean  fluid 
velocity  is  used  in  the  particle  equation  of  motion.  The  influence  of  fluid  turbulence  is  accounted 
for  separately  either  by  adding  a  turbulent  diffusion  force,  ftCj  into  the  particle  equation,  or  by 
correcting  the  particle  velocity  obtained  from  the  single  particle  equation  by  adding  a  turbulent 
diffusion  velocity,  UPd  (Celik,  1988;  Dukowicz,  1980;  Berlemont  et  al.,  1993;  Hallmann  et  al., 
1995).  The  calculation  of  Upa  usually  involves  solving  an  Eulerian  transport  equation  for  the 
particulate  phase.  The  transport  equation  involves  a  transport  coefficient  called  particle 
diffusivity,  which  must  be  prescribed  empirically.  Many  empirical  equations  are  suggested  in  the 
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literature  (e.g.  Mostafa  and  Mongia,  1989)  for  calculating  this  parameter.  The  Eulerian  gas-phase 
equations  contain  source  (or  sink)  terms  to  account  for  the  influence  of  particulate  phase  on  the 
fluid  and  vice  versa.  This  input  can  also  be  computed  from  a  probability  distribution  function 
(PDF)  using  a  Monte  Carlo  method.  The  PDF  is  usually  of  Gaussian  type  with  the  mean  and  the 
variance  being  a  function  of  k  and  s  or  other  turbulence  parameters.  Typically  tens  of  thousands 
of  particle  trajectories  must  be  calculated  to  obtain  reasonable  statistics  for  the  particles.  Particle- 
particle  interaction  is  also  taken  into  account  to  a  certain  extent  using  semi-empirical  relations, 
and  by  allowing  fragmentation,  coagulation  (or  agglomeration). 

RANS  Models 

The  RANS  models  are  derived  by  time-smoothing  and/or  averaging  of  the  local 
instantaneous  equations  which  are  already  volume  averaged  (see  e.g.  Dasgupta  et  al.,  1994, 
Drew,  1992;  and  Morel  and  Bestion,  1997).  Most  researchers  use  some  version  of  the  so  called 
k-e  model  or  algebraic  stress  models  with  marginal  success  (see  for  example  Saif  and  Bertodano, 
1996;  Shimizu  and  Yokomine,  1993;  Sato  et  al.,  1996;  Celik  and  Gel,  1997;  Paterson  et  al., 
1996).  The  central  assumptions  in  the  k-c  model  are  the  same  as  in  single  phase  models,  i.e.  the 
isotropic  eddy  viscosity  concept  and  the  gradient  diffusion  model,  in  addition  to  the  assumption 
of  local  equilibrium  which  is  valid  with  limitations  only  at  very  high  Reynolds  numbers.  As 
reviewed  by  Sato  et  al.  and  Shimizu  and  Yokomine  the  turbulent  kinetic  energy  and  epsilon 
equations  are  most  problematic  for  two-phase  flow  applications.  It  has  been  shown  in  the 
literature  (e.g.  Shimizu  et  al.,  1993;  Balzer  et  al.,  1997;  Celik  and  Gel,  1997)  that  these  equations 
must  be  modified  to  include  the  turbulence  modification  due  to  the  presence  of  dispersed  phase. 
The  void  fraction  fluctuations  can  not  be  neglected  in  the  continuous  phase  equations.  It  is 
particularly  difficult  to  incorporate  higher  order  correlation  terms  which  arise  from  the  particle- 
fluid  interaction  terms,  which  may  account  in  certain  cases  (particle  Reynolds  number  Rep<100) 
for  the  suppression  of  turbulence  and  in  certain  cases  (Rep>  400)  for  enhancement  of  the 
turbulence  of  the  main  stream.  As  demonstrated  by  Celik  and  Gel  (1997)  the  cross-stream 
distribution  of  bubble  concentration  in  shear  layers  is  primarily  determined  by  the  kinetic  energy 
of  the  continuous  phase.  The  more  advanced  Reynolds  stress  models  (see  e.g.  Bertodano  et  al., 
1990;  Simonin  et  al.,  1995)  have  not  proven  to  be  any  more  successful  than  the  k-s  model.  These 
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models  involve  many  more  additional  equations  and  ad-hoc  assumptions  in  the  presence  of  the 
second  phase. 

LES  and  DNS  Approach 

Much  of  the  uncertainties  in  turbulence  modeling  can  be  eliminated  if  the  recently 
developed  Large-eddy  simulation  (LES)  technique  (Galperin  and  Orszag,  1993;  Lesieur  and 
Metais,  1996;  Kamiadakis  et  al.,  1990;  Reynolds,  1989;  Piomelli  et  ah,  1988)  can  be  extended  to 
two-phase  flow  with  appropriate  modifications.  LES  is  a  numerical  technique  in  which  large- 
scale  energy  containing  eddies  (those  responsible  for  the  primary  transport)  are  resolved 
explicitly  and  only  the  small-scale  sub-grid  motions  are  modeled.  This  is  not  the  same  as  the 
Direct  Numerical  Simulations  (DNS)  which  solves  the  full  conservation  equations  without 
parameterization  and  resolves  all  scales  down  to  the  Kolmogorov  scales.  Consequently  it  is 
limited  to  much  lower  Reynolds  numbers  (typically  Reynolds  numbers  less  than  1 04)  and  simple 
geometries  (see  e.g.  Elgobashi,  1996).  For  most  engineering  problems  the  turbulence  Reynolds 
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number  is  in  the  order  great  than  10  (easily  reaching  10  in  geophysical  flows)  which  require 
computations  involving  more  than  10 12  grid  nodes.  This  is  beyond  the  practical  limits  of  present 
computer  technology.  More  over  it  is  also  extremely  difficult  to  extend  the  highly  accurate 
schemes  such  as  spectral  methods  which  are  necessary  for  DNS  to  more  complicated  geometries. 
LES  technique,  on  the  other  hand,  solves  the  unsteady  three-dimensional  Navier-Stokes 
equations  with  an  appropriate  filtering  procedure.  The  filtered  equations  involve  Reynolds  stress 
type  of  terms  (contributions  from  the  sub-grid  scales)  which  are  modeled  by  relating  them  to 
strain  rates  with  an  eddy  viscosity  as  the  proportionality  coefficient.  This  is  referred  to  as  SGS 
closure  model.  The  eddy  viscosity  is  usually  calculated  as  a  function  of  the  mesh  size.  Hence,  the 
finer  the  numerical  mesh  size,  the  less  important  is  the  effect  of  the  smaller  scales  that  are 
filtered  out.  The  three-dimensional  time-dependent  details  of  the  largest  scales  of  motion  (those 
responsible  for  the  primary  transport)  are  computed.  The  size  of  the  scales  that  need  to  be 
resolved  determine  the  numerical  mesh  size  to  be  used.  LES  makes  extensive  use  of 
computational  accuracy  and  computer  power  rather  than  solving  a  large  number  of  empirically 
modeled  equations  as  is  the  case  for  RANS  models.  In  this  regard  it  also  requires  the  use  of 
accurate  numerical  schemes.  In  the  DNS  of  dispersed  two-phase  flows  (see  e.g.  Elgobashi,  1996: 
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Part  1  &  2;  McLaughlin,  1994;  Chen  and  Pereira,  1995;  Rouson  and  Eaton,  1994;  Pedinotti  et  al., 
1993),  in  addition  to  the  Navier-Stokes  equations,  the  instantaneous  particle  equation  is  solved  in 
a  Lagrangian  frame  of  reference  as  explained  above  for  the  EL  approach.  Two-way  coupling  is 
achieved  by  adding  the  momentum  source  terms  in  each  computational  cell  as  the  calculations 
proceed  in  time.  DNS  is  the  most  sophisticated  approach  to  represent  the  details  of  two-phase 
dispersed  flows  at  all  scales.  However,  because  of  the  restrictions  cited  above  the  most  viable 
technique  for  predicting  two-phase  flow  turbulence  is  probably  the  LES. 

LES  has  also  been  applied  (Wang  and  Squires,  1995;  Lavieville  et  al.,  1995;  Ebert  and 
Dehning,  1992)  to  two-phase  flow  turbulence,  but  in  most  studies  only  the  continuous  phase 
turbulence  is  resolved  using  LES,  and  almost  exclusively  the  applications  are  in  the  area  of  gas 
solid  flows  with  the  exception  of  Lapin  and  Lubbert  (1994).  The  later  authors  have  reported 
successful  results  for  gas-slurry  column  reactors.  Their  approach  might  be  classified  as  LES  but 
they  do  not  refer  to  the  LES  literature.  In  gas-solid  flow  applications  the  particles  are  treated 
usually  as  passive  agents  (one-way  coupling)  and  their  motion  is  calculated  using  the  BBO 
equation  in  which  the  filtered  continuous  phase  velocity  components  are  used  directly.  The 
particle  equation  of  motion  is  not  filtered  which  is  inconsistent  with  the  formulation  of  the  LES 
equations.  If  the  filtered  values  are  used  in  the  particle  equation,  the  contribution  from  sub-grid 
turbulence  must  be  included.  The  prediction  of  turbulence  augmentation  by  the  presence  of 
dispersed  phase  (two-way  coupling)  is  not  trivial  and  there  remain  many  issues  to  be  resolved. 

Inflow  Boundary  Conditions 

Although  advantageous  to  RANS  based  turbulence  models  there  are  still  some 
unresolved  issues  in  LES  applications  to  turbulent  flow  simulations.  One  of  the  most  critical 
problems  is  that  large  eddy  simulations  of  spatially  inhomogeneous  flows  require  unsteady 
turbulent  inflow  boundary  conditions.  Especially  if  high  Reynolds  number  flows  exist,  the  inflow 
turbulence  cannot  be  ignored  and  its  effects  has  to  be  accounted  for  via  methods  of  generating 
turbulence  with  prescribed  turbulence  intensity  at  the  inflow  boundary.  Two  main  groups  of 
approaches  in  large  eddy  simulation  studies  are  usually  adopted  to  generate  inflow  turbulence. 
The  first  one  is  to  conduct  the  auxiliary  simulations  of  turbulent  flow  fields  using  the  LES 
approach  (e.g.  Lund  et  al.,  1998  and  Spalart,  1988)  and  to  store  the  time  series  of  fluctuating 
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velocity  components,  which  in  turn  will  be  reintroduced  as  inflow  boundary  conditions  of  the 
main  simulations.  The  main  disadvantage  of  this  approach  is  the  computational  cost  which  has  to 
be  invested.  The  second  group  artificially  generates  time  series  of  random  velocity  fluctuations 
by  performing  an  inverse  Fourier  transform  for  prescribed  spectral  densities.  The  main 
disadvantage  of  this  method  is  that  the  generated  velocity  field  may  lack  both  turbulent  structure 
and  non-linear  energy  transfer.  This  may  be  remedied  by  Celik  et  al.’s  (1999)  approach,  which 
suggest  a  relatively  simple  random  flow  generation  (RFG)  algorithm,  which  can  be  used  to 
prescribe  inlet  conditions  for  spatially  developing  inhomogeneous  anisotropic  turbulent  flows. 

1.3  Objective  and  Present  Contribution 

The  objective  of  this  study  is  to  develop  theoretical  and  numerical  capabilities  and  use 
these  to  study  the  physics  of  bubble  dynamics  and  mixing  in  three  dimensional,  turbulent  flows 
encountered  in  the  near  wake  (~l-3  ship  lengths)  of  ships  using  LES.  Our  focus  will  be  on 
prediction  of  bubble  size  and  concentration  variations,  modeling  of  cluster  formation 
mechanisms  in  vortical  structures,  including  their  influence  on  the  turbulence  of  the  carrier 
liquid.  In  the  first  phase  of  this  study  the  surface  waves  will  be  imposed  at  most  as  boundary 
conditions  a  priori.  The  phenomenon  related  to  wave  breakup,  and  air  entrainment  at  the  free 
surface  by  waves  will  not  be  considered.  The  conditions  at  the  inlet  data  plane  will  be  obtained 
from  other  researchers  who  already  have  grants  from  ONR  to  study  ship  hydrodynamics.  Our 
theoretical  effort  will  be  supplemented  by  laboratory  experiments  which  will  be  used  for  model 
validation.  We  shall  develop  parallel  algorithms  for  high  performance  computers  to  be  able  to 
compute  realistic  particle  statistics  with  the  goal  of  at  least  105  particle  trajectories.  Our 
contribution  will  be  in  advancement  of  the  LES  technique  and  its  application  to  dispersed  two- 
phase  turbulent  flows  with  two-way  coupling.  This  work  will  be  unique  in  bringing  a  new 
formulation  regarding  the  subgrid-scale  turbulence  which  will  be  based  on  a  one-equation  model 
involving  the  kinetic  energy  of  small  scales.  This  equation  will  be  modified  to  account  for  the 
suppression  and/or  enhancement  of  turbulence  by  the  dispersed  phase.  A  new  sub-grid  length 
scale  parameterization  will  also  be  developed.  As  such  this  study  will  be  the  first  (to  the  best  of 
our  knowledge)  attempt  to  apply  the  LES  technique  to  bubbly  ship  wakes. 
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2  METHODOLOGY 


We  performed  this  study  utilizing  a  readily  available  computer  code  obtained  from 
Stanford  University,  with  appropriate  modifications  to  make  the  code  suitable  for  gas-liquid  two 
phase  flows.  The  computer  code  is  based  on  an  essentially  non-staggered  grid,  finite  volume 
method  using  a  fractional  time  step  approach.  Non-orthogonal  curvilinear  coordinates  are  applied 
with  an  overall  second  order  accuracy  in  both  space  and  time.  The  Crank-Nicolson  discretization 
scheme  has  been  applied  for  diagonal  viscous  and  diffusion  terms  and  an  explicit  Adams- 
Bashforth  scheme  is  employed  for  other  terms.  The  central  differencing  (CD)  scheme  (with 
special  care  due  to  numerical  instabilities)  is  applied  to  discretize  the  convective  terms.  Detailed 
information  can  be  found  in  Shi  (2001).  It  has  been  modified  by  our  group  to  run  in  LES  mode 
with  particles  on  parallel  machines  (see  Yavuz  et  al.,  2004;  Shi  et  al.,  2006;  Smirnov  et  al., 
2005). 

The  current  solution  method  is  2nd  order  accurate  both  in  space  and  time.  Theoretically 
the  higher  the  order  of  the  numerical  scheme,  the  better  the  resolution  should  be  under  the  same 
grid  spacing  (Shi,  2001).  However,  Rai  and  Moin  (1991)  have  shown  that  higher  order  of 
accuracy  combined  with  coarse  grid  spacing  does  not  necessarily  give  better  results.  Jordan 
(1999)  showed  that  the  results  could  be  improved  by  improving  the  grid  spacing.  The  power  law 
scheme  is  inaccurate  under  some  limitations,  in  that  when  convection  is  dominant;  it  reduces  to 
1st  order  upwind  scheme  (Patankar,  1980).  On  the  other  hand,  higher  order  CD  schemes  have  in 
addition  the  problem  of  artificial  high  frequency  oscillations  that  may  contaminate  the  turbulence 
field  (Rai  and  Moin,  1991).  In  LES,  explicit  schemes  are  preferable,  but  if  stability  is  an  issue, 
some  implicitness  can  be  introduced,  i.e.  Crank-Nicolson  time  splitting.  For  information  on  the 
LES  code,  i.e.  equations,  time  advancement  and  spatial  discretization  schemes,  the  fractional 
step  method,  poisson  solver,  etc.,  the  reader  is  referred  to  Cehreli  (2004)  and  Shi  (2001). 

Subgrid-scale  closure  for  two  phase  bubbly  turbulent  flows  have  been  considered  at 
various  levels,  starting  from  the  simplest  zero-equation  Smagorinsky  closure,  and  extending  to  1- 
equation,  hybrid  Smagorinsky-turbulent  kinetic  energy  (TKE)  closure,  and,  if  necessary,  to 
dynamic  SGS  eddy  viscosity  model  (for  a  review  see  Ferziger,  (1993)).  The  basic  Smagorinsky 
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model  and  its  variants  need  empirical  modifications  for  wall  damping  and  low  turbulence 
Reynolds  numbers  (see  Reynolds  (1989)).  The  RNG  (Yakhot  and  Orszag,  (1989))  method 
produces  a  modified  form  of  the  Smagorinsky  model  by  an  analytical  scale  elimination  method, 
and  it  leads  to  an  automatic  correction  for  low  turbulence  Reynolds  numbers.  The  need  for  wall 
damping  was  satisfied  using  a  van  Driest  damping  function.  We  conducted  a  comparative  study 
of  various  subgrid-scale  models  and  made  improvements  to  better  model  the  dispersed,  gas- 
liquid  two  phase  flows  relevant  to  ship  wakes. 


2.1  Governing  Equations  and  Navier  Stokes  solver 


The  LES  code  used  was  originally  developed  by  Zang  et  al.  (1994).  The  equations  for  an 
incompressible,  viscous  fluid  flow  in  Cartesian  (physical)  space  can  be  presented  in  terms  of  the 
Cartesian  velocities  Uj  as 
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u 


An  additional  equation  that  represents  the  conservation  of  a  scalar,  such  as  kinetic 
energy,  temperature,  etc.  is, 


d^  d[<puj)  d2  (p 


*  + 

dt  dx : 


=  a 


dx2 


+  Sp) 


(2.3) 


where  a  is  a  material  coefficient  that  could  be  thermal  diffusivity,  conductivity  or 
viscosity,  depending  on  which  scalar  equation  is  solved.  S^tp)  is  a  sink/source  term.  In  the 

above  equations,  Uj  is  the  Cartesian  velocity  vector,  P  is  the  total  pressure,  and  i  ,  j,  k,  are  the 
notations  that  represent  the  directions;  xi  is  the  axial  coordinate,  X2  is  the  vertical  coordinate,  and 
X3  is  the  transverse  coordinate  in  the  Cartesian  coordinate  system.  Qs  is  the  angular  velocity  of 
the  system  rotation  due  to  the  turning  of  the  ship,  v  is  the  kinematic  viscosity  and  by  is  the 
Kronecker  delta  symbol.  In  Eqn.  2.2,  third  term  in  the  right  hand  side  (RHS)  represents  the 
Coriolis  force  and  the  fourth  term  represents  the  Centrifugal  force  due  to  the  turning  of  the  ship. 
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It  should  be  noted  that  Einstein  summation  rule  applies  to  repeated  indices  except  for  the  term, 

The  Navier  Stokes  equations  have  been  developed  for  a  general  purpose  model  of  fluid 
flow  from  basic  principles  of  conservation  of  mass  and  momentum  for  a  Newtonian  fluid. 
Moreover,  the  numerical  models,  such  as  finite  volume  method  (Ferziger  and  Peric,  1997)  also 
incorporate  the  law  of  conservation  of  mass  and  momentum  for  space  integration. 

The  filtering  process  plays  an  important  role  in  distinguishing  small  scales  and  large 
scales  in  LES.  A  flow  variable,  f  can  be  decomposed  into  a  large  scale  of  the  flow  field 
component  that  is  resolved,  /  and  a  small  scale  component  that  is  filtered  out,  /’ ,  as, 


/  =  /  +  /' 


(2.4) 


Then  the  resolved  scale  field  is  obtained  by  applying  spatial  filtering  that  can  be  generally 
expressed  by  the  convolution  integral  (Leonard,  1 974)  for  the  calculation  domain,  D,  as, 

f(xl  x2x3)  =  in  Gj(Xj,x'j  :  A) / (x]x2x3)dx]dx2dxi  (2.5) 

D  7=1 

where  G  is  the  filter  function  and  A  is  the  filter  width,  i.e.  the  wavelength  of  the  smallest 
scale  retained  by  the  filtering  function.  The  most  commonly  used  filter  functions  are  the  box 
filter,  the  sharp  Fourier  cutoff  filter  and  Gaussian  filter  (best  defined  in  wave  space),  and  the  top 
hat  filter  (in  real  space)  (Piomelli,  1999).  In  the  present  finite  volume  formulation,  a  volume 
average  box  filter  used  by  Deardoff  (1970),  is  used,  in  which  Gj=T  (Zang  et  al.,  1993). 

Applying  the  filtering  operator  to  the  governing  equations  and  following  the  formulation 
of  Zang  et  al.  (1993)  in  a  conservative  manner,  the  spatially  filtered  flow  conservation  equations 
can  be  written  as, 


dx, 


=  0 


(2.6) 


^A=s, 

dt  dxj 


(2.7) 
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(2.8) 


c<t>  t  dRj 

dt  8xj 


=s,(*) 


where 


du. 


Flj=uiuJ  +  p5iJ-v^  +  rij 


S,  =  2Q 


fix',  _  fix,  ^ 

j  uix~ui* 

\  dxtJ 


d  -J  <¥ 


(2.9) 


(2.10) 


(2.11) 


Here,  both,  the  SGS  stress,  rlJ  and  the  SGS  flux,  Zj  represent  the  effect  of  the  SGS 
motion.  They  arise  due  to  the  filtering  of  the  nonlinear  advective  terms.  The  formulations  are, 


Zij  ~  UjUj  ~  UjUj 


Zj=uJ<t>-uJ<f> 


(2.12) 


(2.13) 


Hence,  the  SGS  stress  and  the  SGS  flux  both  contain  the  interaction  of  subgrid  scales 
with  themselves  and  with  the  resolved  scales.  In  the  above  equations,  uj  is  the  filtered  velocity 
vector  and  p  is  the  reduced  dynamic  pressure  in  which  the  total  pressure,  P,  is  calculated  as; 

p  =  P«P+x-Paay  (2.14) 

where  r2  is  the  square  of  the  distance  to  the  rotation  axis,  in  terms  of  both  the  axial  and 
the  transverse  directions,  i.e.  r  (=  x{  +  x3 )  and  p0  is  the  reference  density.  It  should  be  noted 

that  the  non-inertial  effects  are  split  with  Coriolis  terms  appearing  as  a  source  term  in  Equation 
2.7  and  since  the  centrifugal  force  term  is  independent  of  the  fluid  motion,  the  effect  of 
centrifugal  force  is  included  in  the  total  pressure  term,  i.e.  Equation  2.14. 
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Since  the  discrete  solution  represents  the  resolved  field  which  is  topped  by  an  overbar 
computed,  thus  the  stress  and  the  flux  terms  should  be  modeled  using  the  resolved  quantities. 
Most  SGS  models  for  ty  and  %  j are  eddy  viscosity  models  of  the  form: 


-2vA  +Cr 


f  8 

Zm _ y_  jm 

ij  0  ^kk 

v  3 


(2.15) 


Xt 


■■-Or^+CP. 


'  J 


(2.16) 


where  v,  is  the  turbulent  eddy  viscosity,  it  represents  the  effects  of  sub  grid  turbulence;  in 
our  case,  a ,  is  the  turbulence  diffusivity  and  SI  is  the  large  scale  (resolved)  strain  rate  tensor, 
defined  as, 


Sj=- 
9  2 


-  \ 


dut  dUj 
KdxJ  + 


dx , 


(2.17) 


'  J 


Eq.  2.15  and  2.16  introduces  two  sets  of  additional  terms  to  the  filtered  governing 
equations:  Cr  and  Lmij  ;  Cr  and  Pj.  Lmij  ,  the  modified  Leonard  term  and  Pj  are  defined  by, 


-  UjUj 

(2.18) 

Pj  =  urf  -urf 

(2.19) 

The  modified  Leonard  term  or  Pj  represent  the  interactions  between  resolved  scales  that 
result  in  sub-grid  scale  contributions  and  can  be  computed  directly  from  the  resolved  flow  field 
(Piomelli,  1999).  The  value  of  the  scale  similarity  coefficient,  Cr  in  Eqn.  2.15  is  either  0  or  1  or 
may  be  determined  dynamically  depending  on  the  type  of  sub-grid  scale  (SGS)  model  being 
used.  When  Cr  =  0,  Eqn.  2.15  represents  the  Smagorinsky  model.  When  Cr  =  1,  it  represents  the 
dynamic  mixed  model  of  Zang  et  al.  (1993). 

To  tackle  problems  of  complex  geometries,  the  above  mentioned  equations  are 
transformed  from  the  physical  to  computational  space  and  formulated  for  a  generalized 
curvilinear  coordinate  system.  The  solution  of  numerical  problems  in  complex  domains  using 
boundary-fitted  curvilinear  coordinates  is  now  a  typical  technique.  The  physical  space  is  denoted 
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(x  x  x  ^  ( E  E  E  ^ 

by  coordinates  '  15  2’  3 ''and  the  computational  space  by  v’1’b2’b3/.  The  chain  rule  of 

derivatives  has  been  applied. 

a  at,  a 


dxj  dxj  84z 


(2.20) 


In  order  to  use  the  finite  volume  discretization,  it  is  desirable  to  cast  the  equations  in  the 
“Strong-Conservation-Law  Form”  as  explained  briefly  in  Zang  et  al.  (1994). 

Substituting  Eqn.  2.15  into  Eqn.  2.7,  and  applying  coordinate  transformation  and 
combining  terms  accordingly,  Eqs.  2.6,  2.7  and  2.8  in  time-dependent  boundary-fitted  curvilinear 
coordinates  are, 


dU, 


K, 


m  _ 


=  o 


d(j  y  )  Qp_ 

_ V _ _j_  im  _  g 

dt  d4m  ' 

8Rm  —  (  .v 
- i  +  —BL  =  SnU) 

8t  84m  TlKf) 


(2.21) 


(2.22) 


(2.23) 


where 


F,m  =  Umul+J-'^p-(v  +  vr)G" 

dxi  d$n 


(2.24) 


St  =2/”1Q5-^l 
dx, 


_  dx 
u 


3  -  Sx, 


V 


+ 


mj 


j- 1  dvr  c  JL 

ax .  ax;.  a^m  a^  r  a^ 


j 


-1 

ax  L,J 

UXj  J 


(2.25) 


Rm-UJ-(a  +  aT)Gn 


8<j> 

Hn 


(2.26) 


.  dxJ  J 


(2.27) 
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cOc 

where  for  the  curvilinear  space  variables;  the  inverse  Jacobian,  defined  as  J~x  =  det  —A  ;  the 
contravariant  velocity  ,  JJ  =  J'1  —^-u  ;  the  contravariant  volume  metrics,  Gmn ,  that  measures 

m  dxj  J 

the  skewness  of  a  grid  cell,  is  defined  as  Gmn  =  .  5'(^),  the  transformation  of 

dxj  dxj 

St  (</>)  to  the  computational  domain,  must  be  changed  accordingly,  depending  on  which  scalar 
equation  is  solved. 

If  the  flux  terms  dFimjd^m  and  dRm/d are  split  like  in  Zang  (1993)  as, 


=  -[Cl+B,{p)  +  DE(ul)  +  Dl(u,j] 


(2.28) 


=  -[e,  +  fe(p)+f,(p)\ 


(2.29) 


where  Cj  and  Ep  represent  the  convective  terms,  Bi  is  the  discrete  operator  for  the  pressure 
gradient  term,  De  and  Di  (Fe  and  Fi)  are  discrete  operators  for  the  explicitly  treated  off  diagonal 
terms  and  the  implicitly  treated  diagonal  viscous  (diffusive)  terms. 


2.2  Numerical  Method 

The  computer  code  is  based  on  an  essentially  non-staggered  grid,  finite  volume  method 
using  a  fractional  time  step  approach.  A  staggered  grid  method  in  curvilinear  coordinates 
requires  a  large  amount  of  computer  memory  for  the  metrics  (Zang  et  ah,  1993),  hence  the  non- 
staggered  method  originally  developed  by  Rhie  and  Chow  (1983)  has  been  used  to  avoid  these 
kind  of  difficulties.  Cartesian  variables  such  as  velocity  and  pressure  are  stored  at  cell  centers 
whereas  the  contravariant  volume  fluxes  are  defined  at  cell  faces  in  a  manner  analogous  to  the 
staggered-mesh  system.  The  volume  fluxes  are  not  solution  variables,  but  rather  are  determined 
through  interpolation  of  the  cell-centered  velocity  values  plus  a  projection  operation  that 
guarantees  exact  conservation  of  mass.  A  traditional  non  staggered  method  does  not  enforce 
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mass  conservation  in  the  cell  and  causes  the  pressure  field  to  decouple  (it  produces  spurious 
oscillations  in  the  pressure  field,  i.e.“checkerboard”  pattern)  (Zang  et  al.  1994),  whereas  the 
method  of  Rhie  and  Chow  (1983)  prevents  the  decoupling  in  its  structure  by  defining  the  volume 
flux  on  its  corresponding  face  of  the  cell  in  addition  to  the  Cartesian  velocity  components  at  the 
cell  center,  therefore  the  momentum  and  continuity  are  both  enforced  in  the  same  control  volume 
and  the  solutions  are  free  from  spurious  pressure  oscillations.  It  is  directly  applicable  to  curved 
domains,  as  the  accuracy  of  the  method  is  not  affected  by  grid  orientations  because  of  the  non- 
staggered  grid  layout.  However,  this  process  eliminates  odd-even  decoupling  at  the  cost  of 
introducing  implicit  4thorder  dissipation,  which  in  turn  may  affect  mass  conservation  (Paterson, 
2003). 

Non-orthogonal  curvilinear  coordinates  are  applied  with  an  overall  second  order  accuracy 
in  both  space  and  time.  The  Crank-Nicolson  discretization  scheme  has  been  applied  for  diagonal 
viscous  (Di)  and  diffusion  (Fi)  terms  in  order  to  remove  the  viscous  instability  (Zang  et  al.,  1994) 
while  an  explicit  Adams-Bashforth  scheme  is  employed  for  all  the  other  terms.  The  off  diagonal 
viscous  terms  (DE)  are  treated  explicitly  in  order  to  simplify  the  LHS  matrix  of  the  momentum 
equation.  The  result,  Zang  (1993),  is 


^4 , 


=  o 


(2.30) 


J 


t(url-»")  3  r  /  .  _  -i  1  r  ,  x  _ 

■' 1 - 1 = |  [cr + de  (k  ) + s-  ]  ~{[crl + de  ( ur 1 ) + sr1 


At 


J 


/  —  n+1  —n  \ 

A  Pi  n 


(2.31) 


At 


2  [s; + fe  [p") + s;  + Fb  (p"'1 ) + s;-'  \ 


(2.32) 


1 

+  — 
2 


where  5/8^  represents  discrete  finite  difference  operators  in  the  computational  space, 

superscripts  represent  the  time  step,  C,  represents  the  convective  terms,  B,  represents  the  dicrete 
operator  for  the  pressure  gradient  terms. 


15 


The  central  differencing  (CD)  scheme  (with  special  care  due  to  numerical  instabilities)  or 
Quadratic  Upstream  Interpolation  for  Convective  Kinematics  (QUICK)  that  calculates  the  face 
value  from  the  nodal  values  using  a  quadratic  upwind  interpolation  is  applied  to  discretize  the 
convective  terms  (Ci).  The  accuracy  of  QUICK  has  been  compared  to  CD,  the  1st  order  upwind 
scheme,  a  hybrid  scheme,  with  the  result  obtained  that  QUICK  produced  good  results  (Zang, 
Hayasa,  1999).  The  spatial  derivates  are  computed  by  2nd  order  central  differences  in  the 
momentum  equations.  Only  the  convective  term  in  the  scalar  equation  (Ep)  is  discretized  using 
the  SHARP  scheme  (Leonard,  1988)  since  it  is  computationally  more  expensive  than  QUICK. 

Because  there  is  no  explicit  equation  to  solve  for  the  pressure  in  time,  the  fractional  step 
method  is  applied  to  solve  the  incompressible  Navier-Stokes  equation.  The  fractional  step 
approach  (Kim  and  Moin,  1985)  or  projection  method,  basically  a  three  step  predictor  corrector 
method,  splits  the  numerical  operators  and  achieves  velocity-pressure  decoupling.  The 
intermediate  velocities  are  interpolated  onto  the  faces  of  the  control  volume  to  form  the  source 
terms  of  the  pressure  Poisson  equation.  The  pressure  field  is  obtained  by  solving  the  pressure 
Poisson  equation  iteratively  with  a  multigrid  method  (Brandt,  1977).  The  true  velocity  field  is 
then  obtained  by  correcting  the  predicted  velocity  with  pressure.  The  steps  are  summarized  from 
Zang  (1993)  as, 

1 .  Predictor  step: 


U-^rAX»;-5i”)  = 


^l|[c;+Ds(i7”)+s"]-i[cr,+fl£(5:"-')+5r']+A(5;")}  e-33) 


2.  Computing  the  pressure  field,  i.e.  finding  <j>; 

If  the  corrector  step  of  the  fractional  step  method  (Equation  2.37)  to  the  Cartesian 
velocity  components  defined  on  a  certain  face  of  the  control  volume, 


u, 


-n+ 1 


U:  -  At 


a* 


(2.34) 
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Combining  Um-Jx  —^-u  with  Equation  2.34,  the  equations  for  £/”+1 , 

dx,  J 


Uf  =U*-At\Gmn^ 

{  K 


(2.35) 


where  U*  =J  ^-u  is  called  the  intermediate  volume  flux.  Since  the  intermediate  velocity 

m  Sxj  J 

u*j  is  defined  at  the  cell  center,  while  the  fluxes  U*m  and  t/”+1  are  defined  on  the  cell  faces,  uf 

% 

has  to  be  interpolated  onto  the  cell  faces  in  order  to  compute  U m . 

By  substituting  Equation  2.35  into  Equation  2.30,  the  pressure  passion  equation  for  <//Hl 
is  obtained  as, 


s  fc_  <y" )  i  { sui 
Si  A  SL) 


(2.36) 


3.  Corrector  step: 


— n+l  *  At  f  n  /  ^n+1  \ 

Ut  ~ui  =-prL5'l^  j. 


(2.37) 


where  /  is  the  identity  matrix,  Uj*  is  the  intermediate  velocity  and  the  scalar  <(>  is  related  to  the 
pressure  p  by 


(2.38) 


Detailed  information  can  be  found  in  Zang  (1993). 


2.3  Sub-grid  Scale  Models 

SGS  stresses  that  are  modeled  (see  Eqn.  2.15)  represent  the  effects  of  the  sub-grid  scale 
motion  on  the  resolved  motion  in  that  they  dissipate  the  resolved  energy  or  backscatter  energy  to 
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the  resolved  eddies.  To  predict  the  flow  dynamics  of  the  wake  behind  a  turning  ship,  the  standard 
Smagorinksy  model  (Smagorinksy,  1963),  based  on  Boussinesq  eddy  viscosity  hypothesis,  has 
been  used.  This  model  is  developed  by  assuming  that  the  small  scales  are  in  equilibrium  so  that 
energy  production  and  dissipation  are  in  balance,  moreover  the  small  scales  dissipate  all  the 
energy  they  receive  from  the  resolved  scales.  This  assumption  is  made  to  simplify  the 
phenomena  and  the  algebraic  model  for  the  eddy  viscosity  is, 

(2.39) 

where  the  filter  length  scale,  A  is  the  volume  average  box  filter  used  by  Deardoff  (1970) 
usually  calculated  as  the  geometric  average  of  mesh  spacings  in  the  Cartesian  directions,  defined 

1/3 

as  A  =  ( Ax{ Ax2 Av3  )  in  finite  volume  formulations  (especially  for  anisotropic  grids)  and  Cs  is 

the  Smagorinsky  constant.  This  constant  is  determined  from  the  isotropic  turbulence  decay.  It  is 
interesting  to  note  that  the  filtering  process  that’s  applied  through  control  volume  approach  may 
use  a  significantly  different  length  scales  in  different  directions  due  to  grid  stretching.  To  see  the 
influence  of  grid  aspect  ratio  on  filtering  process,  a  model  equation  is  used  and  different  filter 
lengths  have  been  applied  (Cehreli,  2004).  The  results  showed  that  the  grid  size  in  each  direction 
influences  mostly  the  degree  of  filtering  in  that  direction;  the  influence  on  the  other  direction  is 
much  lesser. 
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3  THE  RFG  METHOD 


In  a  Reynolds  Averaged  Navier-Stokes  (RANS)  turbulence  modeling  approach 
information  about  turbulent  fluctuations  is  contained  in  the  time  averaged  Reynolds  stresses  of 

the  form  .  These  are  obtained  as  an  outcome  of  a  turbulence  model  that  links  Reynolds 

Stresses  to  mean  flow  quantities  (e.g.  k-<?  model),  or  solves  modeled  transport  equations  for 
each  Reynolds  stress  component  (e.g.  Reynolds  Stress  models).  However,  this  is  not  the  case 
when  the  large  eddy  simulation  (LES)  methodology  is  employed  since  the  goal  here  is  to 
explicitly  resolve  the  turbulent  fluctuations.  In  LES  the  inlet  conditions  can  not  be  derived 
directly  from  experimental  results,  because  of  the  unsteady  and  pseudo-random  nature  of  the 
flow  being  resolved,  unless,  off  course,  the  turbulent  intensity  is  zero  at  the  inlet,  which  is  rarely 
the  case.  This  problem  becomes  more  important  for  spatially  developing  turbulent  flows  where 
for  example  the  boundary  or  shear  layer  thickness  changes  rapidly.  In  such  cases  periodic 
boundary  conditions  can  not  be  specified  like  in  the  case  of  a  fully  developed  channel  flow 
(Ravikanth  and  Pletcher,  2000;  Akselvoll  and  Moin,  1995).  A  similar  situation  exists  when 
prescribing  the  initial  conditions  over  the  whole  calculation  domain.  This  can  be  of  importance 
when  the  turbulent  flow  is  not  steady  in  the  mean  (i.e.  non-stationary  turbulence)  and  the 
transients  of  the  flow  are  to  be  resolved.  Even  for  stationary  turbulent  flows,  if  realistic  initial 
conditions  are  not  prescribed,  the  establishment  of  a  fully  developed  turbulence  takes 
unreasonably  long  execution  time.  For  these  reasons  it  is  necessary  to  initialize  the  flow-field 
with  some  form  of  perturbation  to  provide  the  initial  turbulent  conditions.  It  is  important  that  the 
perturbation  be  spatially  correlated,  as  is  the  case  with  the  real  flow.  For  external  flow  problems 
the  turbulent  flow  field  can  be  initiated  simply  by  appropriately  perturbing  the  inlet  flow-field.  In 
this  case  an  accurate  representation  of  temporal  correlations  of  the  flow-field  can  be  important. 
The  inlet  perturbation  propagates  throughout  the  domain  and  helps  trigger  the  turbulence  that  is 
to  be  captured.  Many  applications  of  LES  begin  with  initializing  the  flow  field  to  that  of  a 
previously  obtained  RANS  solution.  A  higher  resolution  grid  is  then  used  with  an  appropriate 
sub-grid-scale  model.  The  Reynolds  stress  terms  provided  by  the  RANS  solution  can  be  used  to 
construct  spatially  and  temporally  correlated  perturbed  inlet  and  initial  conditions.  In  principle  it 
is  possible  to  predict  turbulence  via  LES  technique  by  starting  from  a  quiescent  flow  or  with  the 
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mean  flow  obtained  from  RANS.  Unfortunately,  it  takes  a  very  long  time  for  a  turbulent  flow  to 
develop  spatially  and  temporally.  This  is  especially  true  in  the  case  of  decaying  turbulence  in  the 
absence  of  strong  turbulence  generating  factors  like  walls.  A  reasonably  accurate  approach  to 
this  problem  is  used  in  modeling  of  boundary  layer  turbulence  (Lund,  1998).  It  consists  in 
applying  a  separate  flow  solver  with  periodic  boundary  conditions  to  construct  the  inlet 
conditions  for  the  LES/DNS  solver.  It  provides  well-formed  inflow  conditions  consistent  with 
the  solution  of  the  Navier-Stocks  equation,  which  makes  it  particularly  suitable  for  DNS. 
However,  its  implementation  may  not  be  straightforward  for  the  problems  without  well  defined 
fully  developed  boundary/shear  layers.  For  some  engineering  problems  it  may  also  be  too 
expensive  in  the  usage  of  computer  resources  and  programming  effort. 

To  remedy  this  problem  the  inlet  and  initial  conditions  are  often  viewed  as  consisting  of  a 
mean  component  and  a  randomly  fluctuating  component  with  the  appropriate  statistics.  Most  of 
the  work  done  in  this  direction  is  based  on  simplified  variants  of  a  spectral  method,  in  which 
Fourier  harmonics  are  generated  with  the  appropriate  statistics  and  assembled  into  a  random 
flow-field.  Realistic  turbulence  spectra  can  be  realized  in  this  way.  In  the  work  of  Lee  et  al. 
(1992)  for  example,  a  very  good  representation  of  turbulence  spectra  was  achieved  by  using 
Fourier  harmonics  with  a  random  phase  shift.  This  is  a  rather  efficient  method  to  generate  the 
inflow  turbulence  with  pre-defined  characteristics.  However,  it  does  not  satisfy  the  continuity  of 
the  flow-field,  which  may  be  important  in  diminishing  the  non-physical  transition  region 
between  the  inlet  flow-field  and  the  solution  provided  by  the  Navier-Stokes  solver  inside  the 
computational  domain. 

A  considerable  amount  of  work  in  random  flow  generation  has  been  performed  in  the 
area  of  particle  dispersion  modeling  using  the  RANS  approach  (Zhou  and  Leschziner,  1991; 
Zhou  and  Leschziner,  1996;  Li  et  al.,  1994).  RANS  modeling  produces  smooth  flow  fields, 
which  do  not  accurately  disperse  particles  that  are  embedded  in  the  flow.  To  correct  this 
turbulent  Reynolds  stresses  are  used  to  generate  temporally  and  spatially  correlated  fluctuations, 
such  that  the  resultant  instantaneous  velocity  can  be  superimposed  on  the  particles  to  induce  a 
realistic  dispersion.  A  number  of  approaches  found  in  the  literature  (Li  et  al.,  1994;  Bechara  et 
al.,  1994;  Fung  et  al.,  1992)  are  based  on  a  variant  of  spectral  method  of  generating  an  isotropic 
continuous  flow-field  proposed  earlier  by  Kraichnan  (1970).  However,  this  flow-field  does  not 
satisfy  the  requirement  of  spatial  inhomogeneity  and  anisotropy  of  turbulent  shear  stresses, 
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which  may  be  important  in  realistic  flows.  The  method  of  Zhou  and  Leschziner  (1991)  complies 
with  the  latter  requirement,  but  the  resultant  flow  field  does  not  satisfy  the  continuity  condition 
and  is  spatially  uncorrelated.  For  homogeneous  isotropic  turbulence,  the  initial  conditions  can 
also  be  constructed  as  described  by  Ferziger  (1983).  The  approach  is  based  on  a  vector  curl 
operation  and  forward/backward  Fourier  transforms,  which  require  a  considerable  computational 
effort.  The  extension  of  this  method  to  anisotropic  inhomogeneous  flows  is  not  trivial.  At  least 
one  study  presents  a  successful  application  of  Kraichnan’s  method  to  anisotropic  flows  (Maxey, 
1987).  The  technique  is  based  on  filtering  and  scaling  operations  applied  to  the  generated 
isotropic  flow-field  to  filter  only  the  vectors  with  the  prescribed  correlations.  Again,  the  filtering 
operation  may  be  expensive  computationally.  The  method  presented  in  this  paper  is  different  in 
that  it  is  based  only  on  scaling  and  simple  coordinate  transformation  operations,  which  are 
efficient  and  fast. 

It  is  the  objective  of  this  study  to  formulate  a  relatively  simple  random  flow  generation 
(RFG)  algorithm,  which  can  be  used  to  prescribe  inlet  conditions  as  well  as  initial  conditions  for 
spatially  developing  inhomogeneous,  anisotropic  turbulent  flows.  In  principle  the  same 
procedure  can  also  be  used  for  initializing  direct  numerical  simulations,  but  the  focus  of  our 
study  is  on  LES,  and  particle  tracking  applications.  The  method  takes  advantage  of  the  previous 
studies  in  the  area  of  particle  dispersion  (Li  et  al.,  1994;  Maxey,  1987).  The  RFG  procedure  is 
developed  on  the  basis  of  the  work  of  Kraichnan  (1970),  and  can  be  used  as  an  efficient  random 
flow-field  generator  in  LES  and  in  particle  tracking  (Shi  et  al.,  2000;  Smirnov  et  al.,  2000).  The 
technique  was  validated  on  the  cases  of  boundary-layer  and  flat-plate  shear  layer  flows  and  is 
further  illustrated  on  the  example  of  bubbly  ship- wake  flow  as  one  of  the  most  challenging  cases 
for  LES  and  particle  dynamics  applications.  Performing  LES  of  ship  wakes  is  particularly 
difficult  given  the  fact  that  the  whole  ship  must  be  modeled  to  capture  a  relatively  thin  3D- 
boundary  layer,  preferably  including  the  viscous  sub-layer.  The  boundary  layer  is  the  source  of 
the  flow  dynamics  that  sets  the  initial  conditions  for  the  wake.  A  simulation  that  includes  the 
whole  ship  and  the  wake  would  require  prohibitively  large  computational  resources.  The  needed 
computational  resources  could  be  substantially  reduced  if  the  appropriate  time-dependent  inlet 
conditions  could  be  constructed  at  the  beginning  of  the  wake,  thus  avoiding  the  need  to  model 
the  entire  ship. 
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3.1  Formulation 


To  generate  a  realistic  flow  field  we  propose  a  modified  version  of  Kraichnan’s  technique 
(Kraichnan,  1970).  The  procedure  we  call  RFG  (Random  Flow  Generation)  combines  the 
advantages  of  the  approaches  mentioned  above  and  is  also  computationally  efficient.  It  involves 
scaling  and  orthogonal  transformation  operations  applied  to  a  continuous  flow-field  generated  as 
a  superposition  of  harmonic  functions.  The  procedure  consists  of  the  following  steps. 

Given  an  anisotropic  velocity  correlation  tensor 

r^UiUj  (3.1) 

of  a  turbulent  flow  field  {ui(xJ,t)}i  j=]  3 ,  find  an  orthogonal  transformation  tensor  a .  that 
would  diagonalize  ry  1 . 

(3.2) 

aikakj=SiJ  (3.3) 

As  a  result  of  this  step  both  ay  and  c„  become  known  functions  of  space.  Coefficients 
c„  =  {C| ,  c3 ,  c3 }  play  the  role  of  turbulent  fluctuating  velocities  («' ,  v' ,  w')  in  the  new  coordinate 
system  produced  by  transformation  tensor  ay .  The  new  algorithm  is  as  follows: 

Step  (i):  Generate  a  transient  flow-field  in  a  three-dimensional  domain  {vi(xJ,t)}l .  ,  3 
using  the  modified  method  of  Kraichnan  (1970) 


v,  (x,  0  =  yJ—£  Ip"  cos(kj  xj + °>J)  +  9"  sin(^”  xj  +  coj )] 


7-L 


■i,  t=~,  c  —  kn  =  k" 

l  T  TJ 


P>£tjmCjK>  (3-6) 

C^>»^(0,1),  K  eA (0,1/2) 

where  l.r  are  the  length  and  time-scales  of  turbulence,  sijk  is  the  permutation  tensor  used 
in  vector  product  operation  (Spain,  1965),  and  N(M,cr)  is  a  normal  distribution  with  mean  M 


1  f.  =  .  Repeated  sub-indexes  imply  summation,  parentheses  around  indexes  preclude  summation. 


22 


and  standard  deviation  a .  Numbers  k" ,  cori  represent  a  sample  of  n  wave-number  vectors  and 
frequencies  of  the  modeled  turbulence  spectrum 

E(k)  =  16(— )1/2ft4  exp(-2&2)  (3.7) 

7t 

Step  (ii):  Apply  a  scaling  and  orthogonal  transformations  to  the  flow-field  v,  generated  in 
the  previous  step  to  obtain  a  new  flow-field  w, 

w;=c0)v0)  (3-8) 

u,=alkwk  (3.9) 

The  procedure  above  takes  as  input  the  correlation  tensor  of  the  original  flow-filed  rr  and 
information  on  length-  and  time-scales  of  turbulence  (/,r).  These  quantities  can  be  obtained 
from  a  steady-state  RANS  computations  or  experimental  data.  The  outcome  of  the  procedure  is  a 
time-dependent  flow- field  u,(Xj,  t)  with  correlation  functions  equal  to  rtJ  and  turbulent 
length/time  scales  equal  to  l,r .  As  will  be  shown  later  this  flow-field  is  also  divergence  free  for 
a  homogeneous  turbulence  and  to  a  high-degree  divergence-free  for  an  inhomogeneous 
turbulence.  By  virtue  of  Eq;  (3.4),  spatial  and  temporal  variations  of  ui  follow  Gaussian 
distribution  with  characteristic  length  and  time-scales  of  l,r .  Sampling  for  a  different  spectrum 
instead  of  Gaussian  can  also  be  used  in  different  problems. 

Equation  (3.8),  provides  the  scaling,  and  (3.9)  -  the  orthogonal  transformation.  Scaling 
factors  c,  obtained  in  step  2  represent  the  scales  of  turbulent  fluctuations  along  each  axis.  They 

do  not  depend  on  time,  whereas  vectors  v,  and  w,  represent  time-dependent  velocity 
fluctuations.  Both  the  scaling  factors  c,  and  the  transformation  tensor  atj  are  computed  in  step  2 

using  efficient  matrix  diagonalization  routines  that  can  be  found  in  standard  libraries  or 
programmed  specifically  for  a  3D  case  to  boost  the  performance.  By  construction,  the  correlation 
tensor  of  the  flow-field  produced  by  Eq.  (3.4)  is  diagonal 

(3-10) 

The  flow-field  wj  produced  after  the  scaling  transformation  (3.8)  is  divergence  free  for  a 
homogeneous  turbulence  and  nearly  divergence  free  for  an  inhomogeneous  turbulence,  as  is 
shown  by  the  following  estimate 
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w  =c  . .  V .  +  C  V.  V  = 

1,1  /,/  /  /  1,1  /  I,/ 


+?;t”cos(-^*;i+Iil^)]=0  (3.11) 

i  i  yv  ^ ]  Cy  /  'T  Cy  /  r 


=>  w, ,  «  0 


(3.12) 


where  we  neglected  all  derivatives  of  c, ,  which  are  slowly  varying  functions  of  x ,  and 
used  the  relation  of  orthogonality  between  £”  and  p",q” 

KP:=Kq:=o 

which  follows  from  the  way  vectors  p" , "  are  constructed  in  (3.6).  For  a  homogeneous 

turbulence  the  approximate  equality  in  (3.11)  will  become  a  strict  equality,  leading  to  a 
divergence  free  flow-field.  In  the  case  of  inhomogeneous  turbulence  the  justification  for 
neglecting  the  derivatives  of  ci  comes  from  the  fact  that  ct  are  derived  from  the  correlation 
tensor  ry,  which  is  produced  by  the  time-averaging  operation  in  (3.1).  As  a  result  of  this 
averaging  rtj  is  independent  of  time  and  may  be  a  slowly  varying  function  of  space  as  compared 
to  the  fluctuating  velocity  w, 


IMHmI  «kJ  -  (3-i3> 

where  |||  denotes  an  appropriate  function  norm.  Relation  (3.13)  justifies  the  first 

approximate  equality  in  (3.1 1),  leading  to  an  approximate  satisfaction  of  continuity. 

The  orthogonal  transformation  (9)  preserves  the  solenoidal  (divergence-free)  property  of 


the  flow-filed 


ui,i  =  auauwJJt  =  5jkwuk  =  wJtj  =  0 


(3.14) 


where  we  used  relation  (3.12)  and  the  rule  of  transformation  of  derivatives:  f  =  a jJj  . 
Thus,  the  generated  flow- filed  ui  is  divergence-free  within  the  accuracy  determined  by  (3.13). 
At  the  same  time  the  new  flow  field  satisfies  the  anisotropy  of  the  original  flow-field  u(x,t) ,  i.e. 


«/«7  aJnwn  = 


aimajnwm  wn  =aimaJncmcnvm  v„  = 


^imajnCmCn^mn  (Cn )  6; 


(3.15) 
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Where  we  used  relations  (3.3),  (3.8)  and  the  last  equality  was  obtained  by  solving  (3.2) 
for  rtJ.  Thus,  the  obtained  flow-field  w,.(x,/)  is  transient,  divergence-free,  inhomogeneous,  and 

anisotropic  with  the  pre-defined  correlation  coefficients. 

Considering  the  flexibility  and  computational  efficiency  of  the  algorithm,  we  should  note 
that  the  random  spectrum  sampling  in  Eq.  (3.6)  can  be  performed  separately  from  the  actual 
assembly  of  the  vectors  in  Eq.  (3.4).  This  leads  to  a  higher  computational  efficiency,  since  the 
spectrum  sampling  can  be  done  outside  of  the  main  time  iteration  loop  of  the  flow  solver  with 
only  the  assembly  of  fluctuating  velocity  components  left  inside  the  time  loop.  This  efficiency 
comes  at  a  price  of  extra  memory  requirements  for  storing  the  random  sample.  The  size  of  this 
sample  is  equal  to  10 -N ,  where  N  is  the  number  of  harmonic  functions  representing  the 
turbulent  spectrum  in  (3.4),  which  is  independent  of  the  actual  grid  size.  So,  for  N  =  1000  and 
double  precision  arithmetics  only  80  KB  of  computer  memory  will  be  needed  to  store  the 
spectrum  for  any  grid  size.  This  offers  a  flexibility  of  shifting  the  priorities  between  the  high- 
accuracy  spectrum  generation  and  speed. 

The  RFG  procedure  can  be  extended  to  include  the  anisotropy  of  turbulence  length-scale. 
In  this  case  instead  of  using  a  scalar  value  for  /  in  (3.5)  one  can  define  a  vector  lr  An 

appropriate  re-scaling  would  be  necessary  to  preserve  the  continuity  of  the  flow-field.  This  can 
be  done  by  introducing  another  scaling  transformation,  similar  to  (3.8),  which  will  guarantee  that 
the  resultant  flow-field  is  divergence-free. 

It  should  be  noted  that  the  turbulent  flow-field  obtained  by  the  RFG  procedure  does  not 
represent  the  solution  of  a  complete  Navier-Stokes  (NS)  equation,  but  rather  of  a  continuity 
equation  only.  This  is  not  a  severe  limitation,  since  the  procedure  is  used  mainly  in  the  context  of 
unsteady  3D  computations,  like  LES,  to  generate  initial  or  inlet  boundary  conditions.  These 
conditions  are  given  on  three-dimensional  subsets  of  a  four-dimensional  computational  domain: 
two  space  and  one  time  dimension  for  the  inlet  flow-filed  and  three  space  and  zero  time 
dimensions  for  the  initial  flow-field.  As  such  these  flow-fields  do  not  have  to  satisfy  the  NS 
equation,  since  for  an  unsteady  LES  this  equation  is  defined  on  a  four-dimensional,  3-space  x  1- 
time  domain.  However,  the  boundary  conditions  should  be  reasonably  continuous,  so  that  the  NS 
solver  will  gracefully  adjust  the  solution  to  the  boundary  conditions  within  the  computational 
domain.  This  is  exactly  what  the  RFG  procedure  is  designed  to  do.  In  addition  to  this  it  provides 
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the  desired  statistical  characteristics  of  turbulence  at  the  boundaries,  like  anisotropy  and 
inhomogeneity,  which  is  of  importance  for  LES  and  unsteady  particle-dynamics  computations. 

Naturally,  this  approach  to  generate  the  inlet/initial  boundary  conditions  is  an 
approximation  and  should  be  used  only  when  the  statistical  features  of  turbulence  at  the 
boundaries  are  of  special  concern  while  the  solution  of  a  full  unsteady  NS  equation  beyond  the 
given  boundaries  is  a  practical  impossibility. 

3.2  Validation  of  RFG  Procedure 

The  first  test  of  the  procedure  was  for  a  homogeneous  isotropic  flow  field.  The  Fourier 

space  was  sampled  with  1000  wave-numbers  selected  according  to  Eq.  (3.4).  Figures  3.1(a), 

\ 

3.1(b)  shows  the  snapshot  of  a  homogeneous  isotropic  velocity  field.  Fig.  3.1(a)  shows  the 
vorticity  field  in  a  cross-section  of  the  computational  domain,  and  in  Fig.  3.1(b)  the  velocity 
distribution  is  presented.  Statistical  post-processing  of  velocity  correlations  was  applied  to  the 

generated  flow-field  in  order  to  verify  that  the  velocity  field  was  isotropic.  For  this  purpose  a 

■> 

turbulent  flow-field  with  the  characteristic  time  scale  of  10' s  was  generated  from  the  Fourier 
spectral  sample  of  1000  wave- vectors  (Eq.  3.4).  The  fluctuating  velocities  were  sampled  at  the 
rate  of  105  FIz  .  Correlations  of  the  fluctuating  velocity  components  were  computed  at  one  point 
in  space  by  averaging  over  time.  Fig.  3.4  shows  the  behavior  of  the  velocity  correlations  as  a 
function  of  the  averaging  time-interval.  The  figure  indicates  convergence  to  the  values 
corresponding  to  Eq.  (3.10). 

The  procedure  was  next  applied  to  a  homogeneous  anisotropic  flow  field  with  c,  selected 

from  typical  boundary  layer  distributions.  This  type  of  anisotropy  leads  to  higher  amplitudes  of 
the  velocity  vectors  in  one  direction  relative  to  the  other  (Fig.  3.1(c)).  The  procedure  was  also 
used  to  generate  the  flow-field  with  anisotropic  length-scales  (see  comments  in  Sec.  3.1).  In  this 
case  the  length-scale  of  fluctuations  was  selected  differently  in  different  spatial  directions.  This 
produced  a  flow-field  where  the  structure  of  the  velocity  fluctuations  seemed  stretched  in  one 
direction  (Fig.  3.1(d)). 

Next  we  applied  the  procedure  to  the  case  of  a  non-homogeneous  anisotropic  boundary 
layer.  Figure  3.1(e)  shows  a  snapshot  of  the  velocity  magnitude  in  the  three-dimensional 
boundary  layer.  An  additional  empirical  factor  related  to  the  boundary-layer  thickness  was 
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introduced  in  this  case  to  better  account  for  the  intermittency  effects.  Figure  3.2  shows  the 
random  signal  produced  by  the  RFG  procedure  sampled  at  different  locations  above  the 
boundary  plane.  As  can  be  seen  from  that  figure  both  anisotropy  and  inhomogeneity  are  evident 
in  the  fluctuating  components.  Experimental  and  direct  numerical  simulation  (DNS)  data  do  exist 
for  this  flow  field,  providing  both  mean  and  fluctuating  velocity  profiles,  as  well  as  turbulent 
correlations.  The  turbulent  boundary  layer  is  two-dimensional  in  the  mean,  though  turbulent 
fluctuations  exist  in  all  three  dimensions,  i.e.  c,  =  u' ,  c2  =  v' ,  and  c3  =  w'  for  the  axial,  vertical, 
and  tangential  directions,  respectively.  In  addition,  the  correlation  involving  the  axial  and 
vertical  velocity  fluctuations  is  significant.  The  Reynolds  stresses  were  obtained  from  (Hinze, 
1975)  where  the  classical  experiments  of  Klebanoff  (1954)  are  summarized. 

A  number  of  realizations  NT  of  the  boundary  layer  was  computed  using  the  turbulence 


time  scale  of  tlurb  =  1 0~3  5 ,  length-scale  of  llrub=  10"3,  and  a  sample  size  of  1000  harmonic 
functions.  The  wave-vectors  for  these  functions  were  taken  from  a  normal  distribution  with  the 
mean  ~  t~]rb .  The  boundary  layer  thickness  (5)  was  allowed  to  grow  according  to  the  following 
empirical  relation: 


£  =  0.16-x 


\  V  J 


=  0.16  -x-Re 


-1/7 


(3.16) 


where  x  is  the  axial  distance  and  U0  is  the  free-stream  velocity,  which  was  set  equal  to 


1.0  m/s .  The  cross  correlation  (mv)  was  normalized  with  the  friction  velocity  UT ,  which  is  itself 


a  function  of  U0 .  The  boundary  layer  thickness  was  randomly  perturbed  with  a  continuous 

function  using  the  same  spectral  sampling  technique  as  for  the  velocity  fluctuations  to  emulate 
intermittency. 

Fig.  3.3  shows  the  vorticity  field  of  the  generated  boundary  layer  flow-field  compared 
with  LES  data  (Speziale,  1998).  As  can  be  seen  from  the  figure,  by  choosing  the  turbulence 
length-scale  correctly,  one  can  achieve  a  good  resemblance  in  the  flow  structure  simulated  with 
this  semi-analytic  approach  and  a  LES  flow-field. 

In  Fig.  3.4,  the  convergence  of  the  random  sampling  process  is  illustrated.  After  about 
500  samples  the  statistics  converge  to  the  prescribed  stationary  correlations. 


27 


To  compare  the  simulation  results  with  the  experimental  data  the  velocity  profile  along  a 
vertical  line  in  the  center  of  the  axial  plane  was  stored  for  each  simulated  realization  of  the  flow- 
field.  The  profiles  of  the  thousand  time  realizations  were  then  used  to  calculate  the  average 
fluctuating  components  in  each  direction,  as  well  as  the  corresponding  cross  correlations.  These 
are  compared  to  the  original  experimental  data  in  Figure  3.5.  As  can  be  seen,  the  experimental 
data  is  well  reproduced. 

The  divergence-free  property  of  the  generated  flow-filed  was  tested  by  computing  the 
divergence  as  a  function  of  turbulence  length-scale  for  three  cases:  isotropic  velocity  field, 
generated  according  to  the  original  Kraichnan  method  (Kraichnan,  1970;  Li  et  al.,  1994), 
anisotropic  velocity  field,  generated  according  to  the  modified  Kraichnan  method,  using  Eqs. 
(3.4)-(3.5)  with  £"  =  £” ,  and  anisotropic  velocity  field  generated  according  to  the  RFG  algorithm 

based  on  Eqs.(3.2)-(3.9).  For  this  test  case  the  anisotropy  of  different  fluctuating  velocity 
components  was  selected  to  be  given  by  a  ratio:  0.1 : 1.4 : 1  for  c,,c2 ,  and  c3  respectively. 


The  divergence  test  was  done  on  a  cubic  grid.  For  each  grid-cell  the  divergence  was 
computed  as  the  sum  of  fluxes  through  cell  faces.  The  Fourier  space  was  sampled  with  1000 
wave-numbers  selected  according  to  Eq.  3.4.  Fig.  3.6  depicts  the  computed  divergence  as  a 
function  of  the  ratio  of  turbulence  length-scale  to  grid  cell  size.  The  result  represents  an  average 
over  1 0,000  realizations  of  the  flow-field.  As  can  be  seen  from  the  figure  in  all  three  cases  the 
continuity  error  decreases  with  the  increase  of  the  turbulence  length-scale.  This  decrease  is 
considerably  slower  for  the  anisotropic  flow-field  generated  according  to  the  original  Kraichnan 
method  compared  with  the  the  cases  of  isotropic  flow  field  and  anisotropic  flow  generated  with 
REG  procedure.  It  should  be  noted  that  the  theoretical  continuity  error  in  the  isotropic  case  is 
zero.  The  discrepancy  between  this  case  and  the  anisotropic  case  computed  with  the  Kraichnan 
method  is  due  to  the  violation  of  continuity  of  that  method  in  the  presence  of  anisotropy.  In 
contrast,  the  flow-field  produced  with  the  new  RFG  procedure  has  practically  the  same  as  for  the 
isotropic  case.  This  means  that  the  anisotropic  flow-field  generated  by  the  RFG  procedure  is 
essentially  divergence  free.  At  the  same  time  this  shows  the  importance  of  scaling  transformation 
for  k'j  in  (3.5)  for  the  fulfillment  of  continuity.  The  upper  divergence  limit  in  the  figure  occurs 


when  the  grid  cell-size  is  comparable  or  greater  than  turbulence  length-scale  / .  It  is  due  to  the 


integration  errors,  which  in  this  limiting  case  can  be  estimated  from  the  relation 


u,ui 
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Another  validation  of  the  method  was  done  on  a  flat-plate  wake  flow,  which  is  a  well 
documented  case  in  the  literature  (Ramaprian  et  al.,  1981).  In  this  case  the  simulation  starts  from 
a  plane  located  behind  the  plate  in  the  wake  region  (Fig.  3.7).  The  inflow  boundary  is  generated 
using  RFG  with  an  input  from  the  experimental  data  (Ramaprian  et  al.,  1981),  including  mean 
velocity  Um,  the  turbulence  intensities  urm , vrms , wrms  and  the  shear  stresses  uv .  In  this  way  we 

can  provide  realistic  inflow  boundary  conditions  including  the  turbulence  characteristics  created 
by  the  plate,  which  is  an  important  factor  for  LES.  The  length  of  the  plate  is  L  =  1.829 m  .  The 
inflow  boundary  is  located  at  19.5 mm  (x/L-  0.01)  measured  from  the  rear  edge  of  the  plate. 
The  computational  domain  size  is  1.0mx0.2mx0.6m  in  x,  y  and  z  direction,  respectively.  The 
corresponding  grid  sizes  are  82x50x50.  Non-uniform  grid  is  used  in  both  x  and  y  directions 
with  stretching  not  exceeding  three  percent  .  The  size  of  the  smallest  cell  is 
0.6mmx0.2mmx\.2mm  while  that  of  the  biggest  cell  is  5Qmmx%mmx\.2mm  in  x,  y,  z 
direction,  respectively.  Neumann  boundary  condition  is  applied  at  the  outflow  boundary. 
Symmetry  boundaries  are  used  in  a  vertical  direction  and  periodic  boundaries  are  used  in  the 
span-wise  direction.  Central  difference  numerical  scheme  and  Smagorinsky  model  are  applied  in 
this  study. 

Fig.  3.8  presents  the  comparison  of  the  turbulent  characteristics  between  the  simulation 
and  the  experimental  results  at  the  inlet  plane.  They  agree  very  well  except  in  the  central  region. 
The  reason  for  the  discrepancies  is  that  the  grid  is  relatively  coarse  in  the  center,  which  leads  to 
the  smoothing  of  very  sharp  gradients.  The  overall  difference  is  under  1  % .  It  should  be  noted 
that  the  agreement  is  so  good  because  the  comparison  is  given  for  the  inlet  plane  where  the  RFG 
procedure  was  designed  to  reproduce  the  turbulence  quantities  exactly. 

Fig.  3.9(a)  shows  the  energy  spectrum  at  the  inflow  boundary.  A  sharp  cut  exists  at  wave 
length  0.01,  which  matches  the  length  scale  we  selected  for  RFG.  Fig.  3.9(b)  and  Fig.  3.9(c) 
show  the  energy  spectrum  at  x  =  0.16  and  x  =  0.53  respectively,  which  were  produced  after  the 
application  of  LES  inside  the  domain.  From  these  figures,  one  can  see  that  a  good  portion  of  the 
inertial  range  is  captured.  With  the  grid  becoming  coarser  in  the  wake,  only  large  wave  lengths 
can  be  resolved. 


2In  this  study,  X  represents  stream-wise,  y  -  vertical  and  z  -  span-wise  directions,  respectively 
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Next  we  investigated  the  turbulence  intensities  downstream  from  the  plate.  From  Fig. 
3.10  to  Fig.  3.13  the  Reynolds  stresses  computed  along  y  -  direction  at  different  x  locations  are 
compared  with  the  experimental  results.  The  agreement  looks  good  although  more  quantitative 
studies  would  be  appropriate.  Near  the  very  beginning  of  the  wake  most  of  the  fluctuations  can 
be  captured.  Later  in  the  flow  field,  the  percentage  of  the  resolved  turbulence  is  becoming  small 
due  to  the  coarser  grids.  It  can  be  noted  that  both  v'  and  uv  at  the  second  x-station  are  higher 
than  those  at  the  first  one.  This  effect  has  been  investigated  and  discussed  by  Nakayama  and  Liu 
(1990)  where  they  attribute  it  to  a  near  wake  phenomenon.  Figure  3.14  shows  the  decay  of  the 
turbulent  kinetic  energy  along  the  center  line  of  the  wake.  Most  part  of  the  turbulent  kinetic 
energy  has  been  captured.  We  anticipate  that  with  finer  grid  sizes  the  turbulence  resolution  will 
be  much  better. 

3.3  Application  of  RFG  procedure 

3.3.1  Boundary  Conditions  for  LES/RANS 

As  an  illustration  of  the  technique  we  used  it  first  in  conjunction  with  the  RANS  and  LES 
methods  to  simulate  turbulent  fluctuations  in  a  ship  wake.  The  high-Reynolds  number  character 
of  ship  wakes  (Re  ~107  -108)  makes  it  rather  time-consuming  to  perform  full-scale  LES  of  these 

flows.  In  this  situation  a  combination  of  a  RANS  and  the  RFG  technique  can  offer  an  efficient 
way  to  restrict  the  LES  runs  to  the  wake  region  only. 

Figure  3.15  shows  a  snapshot  of  an  unsteady  turbulent  flow-field  in  the  inflow-plane 
serving  as  inlet  condition  for  LES  of  a  ship  wake.  Initially  a  steady  state  RANS  solution  is 
obtained  (Figs.  3.15(a),  3.15(c)),  providing  turbulent  stresses  rtj  and  time-scales  coi .  Next  the 

RFG  procedure  is  used  to  generate  continuous  time-dependent  inlet  conditions  with  the  given 
turbulence  characteristics  ( rfJ ,  &>,,  Figs.  3.15(b),  3.15(d)). 

In  another  application  (Fig.  3.16),  namely,  turbulent  flow  around  a  ship-hull  was 
produced  as  a  superposition  of  the  mean  flow  velocity,  computed  with  a  RANS  method  (k-e) 
(Larreteguy,  1999)  and  the  fluctuating  velocity  obtained  with  the  RFG  procedure.  This  flow-field 
can  be  used  to  initialize  the  unsteady  LES  and  for  particle  tracking  applications  (see  Sec.  3.3.2). 
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As  noted  earlier,  the  flow  generated  in  both  cases  may  not  necessarily  represent  realistic 
instantaneous  turbulence  vortex  dynamics  at  the  respective  boundaries,  but  statistically  the  flow- 
field  will  reproduce  the  turbulent  shear  stresses  and  time/length  scales  correctly.  Furthermore, 
the  statistics  of  turbulent  fluctuations  imposed  at  the  space/time  boundaries  may  not  necessarily 
be  preserved  by  the  flow-solver  inside  the  domain.  That  is,  there  may  be  a  transition  region 
between  the  boundary  and  the  inside  of  the  domain  where  the  flow-solver  will  adjust  to  the 
boundary  conditions.  However,  if  the  boundary  conditions  were  generated  using  continuous 
functions  and  realistic  turbulence  quantities  and  spectra,  this  transition  region  will  be  small  and 
the  technique  may  still  present  a  good  alternative  to  solving  the  complete  NS  equation  in  a  larger 
computational  domain. 

3.3.2  Particle  Dynamics  Modeling 

Particle  tracking  in  transient  flows  is  usually  a  time-expensive  computational  procedure. 
In  the  situations  when  the  turbulent  flow  is  computed  using  RANS  models  it  is  possible  to 
compute  particle  dynamics  in  a  steady-state  mean  flow  field  and  add  a  fluctuating  component  to 
particle  velocities.  When  LES  technique  is  used  particles  should  follow  a  time-dependent  flow- 
field  and  the  fluctuating  component  should  still  be  added  to  it  at  smaller  turbulence  scales.  In 
both  cases  the  fluctuating  component  is  derived  from  the  turbulence  intensity  and  length-scales, 
provided  by  the  turbulence  model. 

Bubble-dynamics  in  a  turbulent  flow  is  described  in  Chapter  5  in  detail.  Here,  we  briefly 
present  an  illustrative  example  of  particle-tracking  application  in  a  ship  wake. 

Simulations  of  bubbles  in  ship  wakes  requires  account  of  several  processes,  like  drag,  lift 
and  buoyancy  forces,  bubble  dissolution  in  water,  bubble  interaction  with  the  free-surface 
(including  bubble  disappearance  at  the  surface  and  bubble  generation  due  to  air  entrainment).  In 
some  cases,  because  of  uneven  bubble  distribution  (e.g.  local  clustering),  the  coalescence  and/or 
breakup  of  bubbles  may  be  important.  Because  of  this  non-uniformity,  sharp  gradients  in  bubble 
concentration,  and  low  volume  fraction  of  the  bubbles  Lagrangian  approach  to  model  bubble 
dynamics  is  often  preferred.  Compared  to  the  two-fluids  method  (Elghobashi  et  al.,  1993,  1994, 
1996;  Crowe  et  al.,  1998)  the  Lagrangian  approach  requires  less  empiricism  and  is  more  suitable 
for  parallel  implementation.  We  use  a  particle  dynamics  (PD)  algorithm  based  on  efficient 
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particle  tracking,  population  dynamics  and  a  novel  particle  interaction  techniques  (Smirnov  et 
al.,  2000,  2001). 

To  simulate  bubbles  in  a  ship-wake  we  use  the  combination  of  RFG  and  PD  algorithms. 
Fluid  velocity  at  the  location  of  every  bubble  was  approximated  as  a  sum  of  the  mean  fluid 
velocity  obtained  from  the  RANS  calculations  and  the  fluctuating  part  computed  with  the  RFG 
procedure.  Since  RANS  solution  is  given  only  at  the  Eulerian  grid-node  locations  and  bubbles 
follow  Lagrangian  trajectories,  an  interpolation  is  required  to  approximate  the  mean  velocity  at 
bubble’s  current  location.  No  such  approximation  is  necessary  for  the  fluctuating  part,  since  the 
RFG  procedure  defines  a  flow-field  at  every  point  in  space  and  time.  In  the  simulation  the 
bubbles  were  injected  at  a  single  point  close  to  the  ship  hull  where  the  turbulent  kinetic  energy 
was  near  its  maximum  (Fig.  3.17).  Drag,  buoyancy  and  added-mass  terms  were  included  (Crowe 
et  al.,  1998;  Crowe,  1998).  A  total  of  10000  bubbles  of  100  microns  in  diameter  were 
continuously  injected  into  the  domain.  Two  seconds  of  real-time  were  simulated  for  the  ship- 
length  of  6m  traveling  with  the  speed  of  3  m/s .  The  figure  shows  the  tendency  of  particles  to 
agglomerate  in  dense  groups.  The  characteristic  sizes  of  these  groups  are  in  many  instances 
smaller  than  the  grid-cell  size.  This  reflects  the  very  sub-grid  nature  of  the  RFG  method,  which 
enables  to  capture  finer  details  of  particle  dynamics  than  can  be  resolved  on  an  Eulerian  grid. 

3.3.3  Large-Eddy  Simulations  (LES) 

In  Large  Eddy  Simulations  the  RFG  procedure  can  be  used  to  generate  random  inflow- 
conditions  or  serve  as  a  subgrid-scale  model.  There  is  an  extensive  literature  regarding  LES 
techniques  (Piomelli,  1999).  To  reach  a  state  of  developed  turbulence  in  LES  simulations  require 
a  substantial  computational  time.  Regarding  this  there  are  two  important  problems  in  LES  of  a 
high-Reynolds  number  turbulence  that  can  be  solved  with  the  RFG  procedure:  (1)  assigning 
initial  flow-field  distribution  and  (2)  assigning  turbulent  inflow  conditions.  Conventionally,  the 
first  problem  is  dealt  with  by  increasing  the  transition  phase  of  the  simulation  and  the  second  - 
by  extending  the  size  of  the  computational  domain.  Consequently,  the  remedy  comes  at  a  price 
of  a  longer  execution  time  and  higher  memory  requirements.  By  using  a  RFG  procedure  to 
generate  the  initial  conditions  one  can  cut  down  on  the  execution  time  considerably.  For 
stationary  turbulence  the  approximate  nature  of  the  initial  velocity  distribution  with  respect  to  the 
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solution  of  the  NS-equation  is  of  little  significance,  since  these  discrepancies  are  corrected  in  the 
first  few  iterations  of  the  NS -solver. 

The  problem  of  inlet  conditions  could  be  even  more  important,  since  extending  the 
computational  domain  will  increase  both  execution  time  and  memory  requirements  of  the 
simulation.  In  this  case  RFG  can  provide  reasonable  inflow  conditions  with  the  pre-determined 
anisotropy  properties.  Here  again  we  illustrate  these  advantages  on  the  example  of  a  ship  wake. 
A  complete  LES  simulation  of  the  wake  would  normally  require  simulating  the  unsteady  flow 
around  the  ship  hull  and  in  the  wake  region.  Employing  the  RFG  procedure  for  the  inlet 
conditions,  we  can  restrict  LES  run  to  the  wake  region  only.  In  this  case  the  information  on 
turbulence  levels  and  anisotropy  at  the  inlet  plane,  required  by  RFG,  can  be  obtained  from 
relatively  inexpensive  RANS  calculations. 

After  validating  this  approach  on  the  case  of  a  flat-plate  (see  Sec. 3. 2),  we  applied  it  to  the 
wake  of  a  model  ship  (Navy  5415  model  (Carrica  et  al.,  1998)).  As  in  the  flat-plate  case  the 
inflow  boundary  was  constructed  using  the  data  from  RANS  calculations  (Larreteguy,  1999) 
(Fig.  3.15).  The  turbulent  normal  stresses  are  based  on  the  kinetic  energy.  The  time  scale  and 
length  scale  which  are  used  in  generation  of  the  perturbation  at  the  inlet  plane  are  selected  from 
the  corresponding  relation  between  the  turbulent  kinetic  energy  and  its  dissipation  rate  provided 
by  the  RANS  calculations. 

Figure  3.18  shows  the  instantaneous  stream- wise  velocity  contours  and  vertical  vorticity 
contours,  respectively,  in  a  plane  parallel  to  the  free  surface,  where  some  of  the  turbulence 
structures  in  the  wake  can  be  observed.  Small-scale  turbulent  structures  can  be  seen  in  both 
figures  in  the  near  wake  region.  These  structures  tend  to  increase  in  the  far  wake.  This  can  be  due 
to  the  following  two  factors:  (1)  In  the  very  near  wake,  fine  grids  are  applied  so  that  smaller 
turbulence  structures  are  captured.  (2)  Physically,  larger  turbulence  structures  include  more 
energy  so  that  they  can  last  longer,  while  smaller  turbulence  structures  have  less  energy  and  die 
quickly  by  dissipation.  Another  phenomenon  is  the  increase  of  the  width  of  the  wake  in  the 
downstream  direction.  The  mean  velocity  profile  (not  shown  here)  also  supports  this  result. 
Capturing  of  these  phenomena  testifies  to  the  validity  of  our  approach  of  computing  the  inflow 
boundary. 
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3.4  Discussions 


The  analytical  method  of  Kraichnan  (1970)  was  modified  to  account  for  the  effects  of 
inhomogeneity  and  anisotropy  of  turbulent  shear  stresses.  The  technique  was  realized  in  an 
efficient  random  field  generation  (RFG)  algorithm  which  was  tested  on  the  cases  of  isotropic 
turbulence,  channel  flow  anisotropic  flows,  boundary  layer  inhomogeneous  anisotropic 
turbulence  (Celik  et  al.,  1999)  and  flat-plate  wake  flow.  The  simulated  flow  fields  are  transient  in 
nature  and  satisfy  the  conditions  of  continuity  and  anisotropy  for  homogeneous  flows  and  to  a 
good  approximation  satisfy  these  conditions  for  inhomogeneous  flows. 

The  RFG  procedure  offers  a  relatively  inexpensive  way  to  generate  random  velocity 
fluctuations,  representing  a  turbulent  flow-field.  Since  the  generated  velocity  field  satisfies  the 
relations  of  continuity  and  anisotropy  it  is  a  far  more  realistic  representation  of  turbulence  than 
can  be  obtained  with  a  simple  Gaussian  velocity  distribution  using  a  random-number  generator. 
Because  the  flow-field  produced  by  RFG  may  not  satisfy  the  momentum  equations  it  is  still  an 
approximation.  However,  in  some  applications  this  approach  may  offer  a  simple  and  reasonably 
accurate  way  to  model  turbulence  without  solving  the  complete  Navier-Stokes  equation,  which 
would  require  much  more  memory  and  execution  time. 

In  practical  applications  the  RFG  procedure  provides  the  flexibility  of  a  trade-off  between 
the  accuracy  of  representing  a  turbulent  spectrum  and  memory/time  requirement.  By  increasing 
the  spectral  sample  size  N  in  Equation  3.4  one  can  increase  the  accuracy  of  reproducing  the 
turbulent  spectrum  at  the  cost  of  longer  execution  time  and  higher  memory  utilization.  In 
addition  to  that,  since  the  velocity  field  is  calculated  by  analytical  functions  it  is  given  at  any 
point  in  space  and  time,  and  not  just  at  the  grid  nodes  and  at  discrete  time  values.  Because  of  this 
quality,  the  method  has  a  potential  as  a  subgrid-scale  model  for  LES  or  RANS  simulations  and  in 
modeling  turbulent  particle-laden  flows,  although  its  validity  in  this  respect  would  require  a 
separate  study. 

Advancing  the  realistic  LES  run  to  the  stage  of  developed  turbulence  may  require  days  of 
computation  time3.  Similarly,  to  obtain  realistic  turbulent  inflow  conditions  may  require  the 
extension  of  the  computational  domain  with  the  corresponding  increase  in  computer  time  and 
memory  requirements.  The  RFG  technique  can  reduce  the  flow  initialization  time  to  several 


3Our  benchmarking  was  performed  on  the  533  MHz  DEC-Alpha  processor 
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hours  and  can  be  used  to  continuously  supply  the  turbulent  inlet  conditions  close  to  the  domain 
of  interest,  thereby  reducing  time  and  memory  requirements  of  the  LES  simulations. 

This  study  showed  the  feasibility  of  applying  a  hybrid  LES  technique  in  combination 
with  RFG  algorithm  to  high-Reynolds  number  flows,  like  those  of  ship  wakes.  It  is  also  shown 
that  the  technique  can  be  used  effectively  in  conjunction  with  a  Lagrangian  particle  dynamics 
approach,  is  appropriate  for  bubble  tracking  in  the  wake  and  can  be  easily  incorporated  into  LES 
codes. 
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(e)  Fluctuating  velocity  in  the  boundary  layer 
Figure  3.1:  Simulated  flow-field  using  RFG 


(c)  Anisotropic  velocity  (d)  Anisotropic  length-scale 
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(a)  Axial,  y *  —  0.13 


(b)  Vertical,  y*  =  0.13 


(c)  Tangential,  y*  =  0.13 


(d)  Axial,  y*  —  0.46 


(e)  Vertical,  y*  =  0.46 


(f)  Tangential,  y*  =  0.46 


(g)  Axial,  y*  —  0.76  (h)  Vertical,  y*  —  0.76  (i)  Tangential,  y*  =  0.76 


Figure  3.2:  Instantaneous  velocity  vs  time  step  at  different  locations,  y*  =  y/5 


(c)  RFG  (small  length-scale) 


Figure  3.3:  Vorticity  contours  in  the  boundary  layer 
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(b)  Cross  correlations 

Figure  3.4:  Convergence  of  velocity  correlations 
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Figure  3.8:  Turbulence  i 


(b)  Span-wise,  vrms 
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at  the  inflow  boundary. 


Figure  3.9:  Energy  spectrum  at  different  x  locations 
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Figure  3.10.  Comparison  between  simulation  and  measured  turbulence  intensity  ( urms ). 

- ,  present  simulation  results  at  x  =  3 1 .75 mm  ; 

- ,  present  simulation  results  at  x  =  1 58.75mm  ; 

- ,  present  simulation  results  at  x  =  361 .95mm  ; 

. ,  present  simulation  results  at  x  =  590.55mm  ; 

□,  Exp.  results  (Ramaprian,  et  al.  1981)  at  x  =  31.75mm ; 

° ,  Exp.  results  (Ramaprian,  et  al.  1981)  at  x  =  158.75mm  ; 

A ,  Exp.  results  (Ramaprian,  et  al.  1981)  at  x  =  361 ,95mm  ; 

#,  Exp.  results  (Ramaprian,  et  al.  1981)  at  x  =  590.55mm  ; 
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Figure  3.11.  Comparison  between  simulation  and  measured  turbulence  intensity  (vrm ). 


Figure  3.12:  Comparison  between  simulation  and  measured  turbulence  intensity  ( wrms ). 


Figure  3.13:  Comparison  between  simulation  and  measured  shear  stresses 


45 


Figure  3.14:  Kinetic  energy  profile  along  the  center  line  in  the  wake  of  a  flat  plate 
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0 


(b)  Streamwise  (RANS+RFG) 


(d)  Vertical  (RANS+RFG) 

:  inlet  conditions  for  LES  of  a  ship-wake. 


Figure  3.16:  Turbulent  velocity  around  a  ship  hull  computed  with  the  RFG  algorithm.  (View  from 

below.) 


*  *■» 


Figure  3.17:  Bubbles  in  a  ship-wake.  (Background  shadiiig  is  according  to  the  turbulent  kinetic 

energy.) 
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4  SHIP  WAKE  SIMULATIONS 


The  subject  of  this  Chapter  is  the  application  of  the  numerical  schemes  and  SGS  models 
described  in  previous  chapters  to  the  calculation  of  turbulent  ship  wake  flows.  Starting  from  the 
review  of  the  current  research  status  of  ship  wake  hydrodynamics,  this  Chapter  leads  to  the 
summary  of  the  strategies  described  in  previous  chapters.  Several  cases  have  been  simulated  and 
results  are  analyzed  and  assessed. 

Ship  hydrodynamics  has  been  the  subject  of  numerous  studies  in  the  past.  Both 
experimental  and  numerical  achievements  have  been  made.  Since  real  ships  can  not  be  fitted  into 
tow  tanks,  it  is  common  practice  to  base  the  research  on  small-scale  models  of  real  ships. 
Experimentally,  in  the  60’ s  or  70’ s  of  the  last  century,  most  of  the  experiments  concerned  the 
ship  surface  drag  force,  or  propeller  influence.  Recently,  experiments  are  more  focused  on  the 
investigation  of  micro  dynamics  of  flow  around  ships,  such  as  wave  dynamics,  turbulence 
parameters,  bubble  dynamics,  etc.  Although  experiments  can  provide  useful  results,  there  are 
many  limitations  encountered  in  practice.  First,  it  takes  a  very  long  time  to  design,  set,  and 
perform  the  experiment.  Secondly,  it  is  very  expensive.  Thirdly,  it  can  not  provide  very  detailed 
information  due  to  the  limited  space  to  arrange  the  measuring  devices  or  detectors.  Benefited 
from  the  rapidly  developing  capabilities  of  computers,  the  modem  research  in  computational 
hydrodynamics  has  replaced  the  more  elegant,  but  less  general,  analytical  methods  (Larsson  et 
al.,  1998).  Since  ship  wake  flows  are  the  main  subject  of  this  study  a  separate  but  concise  review 
is  presented  next. 

4.1  Experiments  on  Turbulent  Ship  Wake  Flows 

In  a  earlier  study  of  ship  wake  flows,  Naudascher  (1965)  found  that  the  wake  width  has  a 
power  law  behavior.  This  was  supported  by  Milgram  et  al.  (1993).  By  carrying  out  field 
measurements  for  ship  wakes,  they  found  that  the  wake  width  has  a  power  law  behavior  of  the 
type  x1/5  where  x  is  a  distance  from  the  ship  stem.  Reed  et  al.  (1990)  summarized  the 
hydrodynamics  research  of  ship  wake  flows  both  in  experimental  and  numerical  areas  and 
pointed  out  the  hydrodynamic  wake  schematic  shown  in  Fig.  4.1 .  Benilov  et  al.  (2000)  found  that 
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even  in  the  far  wake  the  ship  wake  turbulence  is  still  well  detectable  and  the  Kolmogorov  range 
can  be  identified.  Hoekstra  (1991)  gave  relatively  comprehensive  tests  of  turbulence  parameters 
in  which  the  wakes  of  different  ship  models  were  studied.  The  vortices  generated  by  the  ship 
body,  named  bilge  vortices,  were  analyzed,  and  it  was  found  that  these  vortices  can  be  avoided 
by  hull  form  adjustment  but  this  was  found  to  be  impractical. 

Stem  et  al.  (2000)  present  some  results  for  CFD  validation,  but  the  experiments  providing 
turbulence  intensities  were  still  underway.  For  this  study  it  was  not  possible  to  find  the  necessary 
experimental  information  to  initialize  an  LES  study,  or  to  build  the  inflow  boundary  condition 
including  the  turbulence  features  induced  by  the  ship  hull.  Thus  the  current  simulation  is  based 
on  RANS  results. 

4.2  Computations  of  Turbulent  Ship  Wake  Flows 

/ 

Most  of  the  computational  fluid  dynamics"  efforts  applied  to  flow  past  ships  are  based  on 
Reynolds-Averaged  Navier-Stokes  (RANS)  equations  utilizing  various  turbulence  models 
(Sotiropoulos  and  Patel,  1995;  Paterson  et  al.,  1996;  Ratcliffe,  1998).  The  commonly  used 
models  include  k-e,  k-co  and  algebraic  stress  models.  RANS  is  often  quite  adequate  for  mean 
flow  predictions,  but  provides  only  limited  information  about  turbulence  characteristics  and 
almost  no  details  on  large-scale  unsteady  structures  of  the  flow  field.  The  LES  technique,  on  the 
other  hand,  was  designed  to  capture  the  unsteady  behavior  of  at  least  the  large  coherent  turbulent 
eddies.  There  is  hardly  any  study  in  the  literature  where  LES  has  been  applied  to  flow  around 
ship-hulls  including  the  wake.  The  main  reason  obviously  lies  in  the  fact  that  large  computer 
resources  are  required.  In  LES,  the  energy  containing  eddies  of  the  flow  are  computed  explicitly, 
while  only  the  more  universal  (isotropic)  small  eddies  are  modeled.  Thus  very  fine  grids  have  to 
be  applied  in  order  to  resolve  the  boundary  layer  near  the  wall  where  the  turbulence  length  scales 
tend  to  zero  as  the  wall  is  approached. 

However,  in  some  applications,  like  bubble  dynamics  modeling,  it  is  still  necessary  to 
resolve  coherent  flow  structures  -  large  turbulent  vortices  and  eddies.  In  RANS  simulations, 
especially  those  using  two-equation  models  such  as  k-s  model,  these  unsteady  flow  features  are 
usually  smeared  out. 
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Since  it  can  be  computationally  prohibitive  to  include  both  the  ship  hull  and  the  wake  in 
LES,  it  would  be  desirable  for  the  purposes  of  pure  wake  simulations  to  start  the  computations 
somewhere  in  the  near  wake  excluding  the  ship  hull.  However,  this  technique,  also  known  as  the 
Initial  Data  Plane  (IDP)  approach  (Hyman,  1998;  Paterson  et  ah,  1996;  Dommermuth  et  ah, 
1996)  can  introduce  considerable  errors  (Hyman,  1998)  -  not  surprisingly  though,  since  the  body 
generating  the  wake  is  not  included  in  the  non-steady  LES  simulations. 

In  principle  it  is  possible  to  predict  turbulence  via  the  LES  technique  by  starting  from  a 
quiescent  flow  or  with  the  mean  flow  obtained  from  RANS.  Unfortunately,  it  takes  a  very  long 
time  for  a  turbulent  flow  to  develop  spatially  and  temporally  without  any  initial  perturbation. 

This  is  especially  true  in  the  case  of  decaying  turbulence  in  the  absence  of  strong 
turbulence  generating  factors  like  walls.  A  reasonably  accurate  approach  to  this  problem  was 
used  in  modeling  of  boundary  layer  turbulence  (Lund  et  al.,  1998).  It  consisted  of  applying  a 
separate  flow  solver  with  periodic  boundary  conditions  to  construct  the  inlet  conditions  for  the 
LES/DNS  solver.  It  provides  fully  turbulent  inflow  conditions  consistent  with  the  solution  of  the 
Navier-Stokes  equation,  which  makes  it  particularly  suitable  for  DNS.  This  method  was  used 
later  by  Wang  and  Moin  (2000)  to  compute  the  inflow  boundary  for  simulating  a  trailing-edge 
flow.  However,  its  implementation  may  not  be  straightforward  for  the  problems  without  well 
defined  fully  developed  boundary/shear  layers.  For  some  engineering  problems  it  may  also  be 
too  expensive  in  terms  of  computer  resources  and  programming  effort. 

Another  numerical  simulation  with  time  dependent  turbulent  inflow  boundaries  was 
performed  by  Voke  and  Potamitis  (1994)  for  a  wake  of  a  flat  plate.  In  their  simulation,  the  inflow 
boundary  information  was  derived  from  a  separate  simulation,  called  the  precursor  simulation,  in 
which  a  low  Reynolds  number  turbulent  boundary  layer  was  simulated.  At  the  inflow  plane  of 
the  wake  simulation  the  three  velocity  components  were  specified  by  reading  one  transverse 
plane  of  velocity  data  from  the  boundary  layer  simulations  at  each  time  step.  The  transverse 
plane  was  selected  at  a  location  where  the  flow  was  fully  turbulent  but  with  a  safe  distance  from 
the  outflow  boundary  of  the  precursor  simulation.  Two  separate  sets  of  velocity  data  of  the 
boundary  layer  were  needed  to  provide  inflow  conditions  for  the  upper  and  lower  halves  of  the 
wake  simulation.  To  achieve  this  they  used  the  data  of  the  same  precursor  simulation  of  the 
boundary  layer  at  well-separated  times  in  it’s  history.  This  method  is  costly  and  only  applicable 
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to  parabolic  wake  flows.  In  application  to  the  ship  wake  flow  it  would  require  LES  of  the  flow 
around  the  hull  itself  which  is  impractical. 

To  achieve  a  higher  accuracy  of  ship-wake  simulations  and  still  remain  within  realistic 
constraints  on  computer  resources,  the  whole  problem  of  computing  the  flow  past  ships  and  their 
wakes  is  divided  into  two  parts:  steady-state  RANS  calculations  around  the  ship  hull,  and  non¬ 
steady  LES  of  the  wake.  In  this  approach  the  LES  starts  from  a  plane  in  the  near  wake  where  the 
inflow  conditions  are  retrieved  from  the  RANS  results.  This  is  similar  to  the  IDP  approach  used 
by  Elyman  (1995)  and  Dommermuth  et  al.  (1993)  where  the  data  reflects  the  mean  flow,  mean 
turbulence  quantities  and  the  scalar  field  at  this  plane.  In  the  present  approach,  a  wider  range  of 
turbulent  quantities  at  the  inlet  plane  are  allowed,  which  include  length-  and  time-scales,  and  all 
nine  components  of  the  turbulent  shear-stress  tensor  (only  six  components  to  be  computed  due  to 
symmetry),  thus  providing  anisotropic  and  inhomogeneous  turbulent  inlet  conditions.  The 
technique  is  based  on  a  Random  Flow  Generation  (RFG)  algorithm  (Celik  et  al.,  1999;  Smirnov 
et  al.,  2001b;  Shi  et  al.,  2000b),  which  is  applied  to  generate  the  inflow  turbulence  on  the  basis  of 
the  information  obtained  from  experiments  or  RANS  calculations.  The  generated  velocity 
fluctuations  satisfy  instantaneous  continuity  equation,  and  the  turbulence  statistics  (Reynolds 
stresses)  are  prescribed  a  priori.  Thus  in  some  sense,  although  the  LES  starts  at  a  plane  behind 
the  body,  the  influence  of  the  body  is  implicitly  included.  The  features  of  the  generated  flow- 
field  such  as  continuity,  anisotropy  and  inhomogeneity  make  the  RFG  method  also  well  suited 
for  setting  the  initial  conditions  for  LES. 

In  the  present  simulations  anisotropic,  inhomogeneous,  unsteady  IDP  conditions  were 
applied.  Thus  the  flow  develops  according  to  the  dynamics  prescribed  by  the  Navier-Stokes 
equations  without  forcing  and  it  is  believed  more  reasonable  in  the  near  wake.  This  approach 
differs  from  that  of  Dommermuth  et  al.  (1996)  which  is  appropriate  for  the  far  wake  field 
(Hyman,  2001). 

To  describe  the  dynamic  behavior  of  the  turbulence  in  the  ship  wake,  in  another  study 
(Benilov  et  al.,  2000)  made  the  following  assumptions: 

The  wake  turbulent  kinetic  energy  significantly  exceeds  the  upper  layer  turbulence  that 
reduces  the  turbulent  wake  problem  to  the  turbulent  region  development  in  a  non-turbulent 
liquid. 
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The  main  source  of  turbulence  is  a  moving  ship  which  means  that  all  interactions 
between  the  wake  turbulence  and  environment  do  not  contribute  in  the  wake  dynamics  and  allow 
the  problem  to  be  reduced  to  the  shear-free  turbulence  model. 

Through  all  this  study,  those  assumptions  were  taken  to  simplify  the  problem  which  is  to 
perform  LES  on  the  wake  flow  of  the  ship  model  DTMB  5415  (5512). 

4.3  Ship  Wake  Simulations  on  a  Straight  Track 

The  ship  wake  simulations  in  the  present  study  are  for  the  ship  model  DTMB  5512,  a 
smaller  ship  model  scaled  from  DTMB  5415  which  are  well  known  in  the  ship  hydrodynamics 
community.  Model  5415  is  a  towing  tank  model  representing  a  modem  naval  combatant.  It  was 
constructed  out  of  wood  and  measured  5.72  m  (18.767  ft)  between  perpendiculars.  Model  5512 
measures  3.04  m  (10.000  feet)  between  perpendiculars.  It  is  an  appended  version  of  5415,  with 
shafts,  stmts,  rudders  and  propellers.  However  in  this  study  it  has  been  simplified  as  a  bare  hull 
model.  The  geometry  of  this  ship  model  is  shown  in  Fig.  4.2.  Note  that  this  picture  is  not  scaled, 
and  it  is  only  used  to  show  the  profile  of  the  ship  model. 

The  Reynolds  number  is  4.65x1 06  and  the  Froude  number  is  0.28.  The  length,  L,  of  the 
ship  model  is  3.048m.  The  ship  is  located  on  the  region  of  x/L  =  0  to  x/L  =  1.0.  The  inflow 
boundary  is  located  at  x/L  =  1.05,  a  small  distance  downstream  of  the  rear  end  of  the  ship.  The 
computational  domain  size  is  1.5m  x  0.3m  x  0.6m  in  x,  y  and  z  direction,  respectively,  with  the 
corresponding  grid  size  of  162  x  50  x  66  (Fig.  4.2).  For  convenience,  the  streamwise  direction 
(x)  of  the  computational  domain  is  marked  from  0  to  1.5  (x  =  0  corresponding  to  x/L  =  1.05, 
where  the  inlet  plane  is  located).  Non-uniform  grid  is  used  in  both  x  and  y  directions  with 
stretching  not  exceeding  one  percent.  The  size  of  the  smallest  cell  is  0.0015x0.0028x0.0068 
while  that  of  the  largest  cell  is  0.028x0.007x0.012  in  x,  y,  z  direction,  respectively. 

We  have  also  performed  simulations  with  various  other  coarser  grids.  Results  from  those 
showed  consistent  trends  in  increased  resolution  of  finer  structures  in  the  flow  field  and  turbulent 
kinetic  energy  as  the  grid  cell  size  decreased.  Here  we  present  the  results  from  the  finest  grid 
resolution  we  were  able  to  achieve  on  a  single  processor. 
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4.3.1  Boundary  Conditions 

The  inflow  boundary  is  generated  using  the  RFG  method  in  conjunction  with  the  RANS 
calculations  (Stem  &  Wilson  2000),  providing  the  mean  velocity  Um,  the  turbulence  kinetic 
energy  k  and  the  specific  dissipation  co.  The  time  scale  x  can  be  taken  proportional  to  l/co.  The 
length  scale  can  be  calculated  from  k  and  co.  In  this  study,  the  length  scale  and  time  scale  used  in 
RFG  were  selected  as  constant.  The  length  scale  was  0.02  of  the  ship  length,  and  the  time  scale 
was  0.001,  non-dimensionalized  by  free  stream  velocity  and  ship  length.  These  numbers  were 
selected  so  that  the  turbulence  length  scale  is  about  1 5%  of  the  ship  width.  Neumann  boundary 
condition  was  applied  at  the  outflow  boundary.  Symmetry  conditions  are  used  in  y  direction  and 
periodic  boundary  conditions  are  used  in  the  spanwise  (z)  direction.  At  the  free  surface  slip  in  x 
and  z  directions  is  allowed  but  the  velocity  component  normal  to  the  free  surface  is  set  to  zero. 
As  such  the  free  surface  is  approximated  as  a  moving  flat  plane. 

The  periodic  boundary  conditions  in  the  span-wise  direction  were  selected  for  the  reason 
of  preserving  the  possible  lateral  fluctuations  caused  by  unsteady  motion  of  large  eddies.  It  is 
also  more  exact  than  the  zero  gradient  boundary  condition  if  applied  to  the  boundaries  in 
spanwise  direction.  Any  other  boundary  condition  would  probably  dampen  that  motion  and 
introduce  artificially  high  subgrid-turbulence  levels  at  the  boundaries.  The  boundary  conditions 
at  the  outlet  were  actually  of  a  free-gradient  (zerostress)  type. 

The  numerical  schemes  and  SGS  model  used  in  the  ship  wake  simulation  are  the  same  as 
in  the  flat  plate  wake  (Shi,  2001 ;  Smirnov  et  al.,  2001a). 

4.3.2  Results  and  Discussion 

The  streamwise  vorticity  contours  at  different  planes  are  shown  in  Fig.  4.3  and  the 
contours  of  the  vertical  component  of  vorticity  are  shown  in  Fig.  4.4.  Fig.  4.3  shows  that 
concentrated  vorticity  decreases  with  axial  distance,  and  in  the  far  wake  vorticity  is  only 
concentrated  near  the  free  surface,  while  the  size  of  the  wake  increases  in  axial  direction  of  the 
flow.  It  is  interesting  to  note  the  two  distinctly  concentrated  vorticity  streaks  away  from  the 
center  line  of  the  wake  (Fig.  4.4).  These  resemble  the  bilge  vortexes  (Reed  et  al  1990)  generated 
by  the  ship-hull  Fig.  4.1.  The  streamwise  velocity  contours  at  a  plane  near  the  free  surface  are 
shown  in  Fig.  4.5.  It  is  interesting  to  see  that  the  vorticity  wake  is  slightly  wider  than  that  of  the 
velocity  wake. 
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Figure  4.1:  a)  Hydrodynamic  wake  schematic  showing  hiajor  large-scale  contributors  to  the  wake 


(Reed  et  al.  1990),  b)  The  schematic  of  the  ship-wake 

A  comparison  with  experimental  results  reported  by  Reed  et  al.  for  the  same  ship  model 
(1990)  revealed  a  qualitative  agreement.  In  particular  the  axial  velocity  contours  at  the  cross- 
section  of  x/L  =  0.61  shown  in  Fig.  4.6  have  the  same  structure  as  the  experimental  flow-fields 
presented  in  that  study  (Figs.  36,  37).  The  same  conclusion  C&n  be  drawn  for  the  velocity  vectors 
at  this  cross-section  (Fig.  4.7).  The  resolved  component  Of  turbulent  intensity,  which  in  our 
simulations  is  about  3.5%,  as  can  be  deduced  from  (Fig.  4.12).  agrees  closely  with  the  measured 
total  intensity  at  X/L  =  0.61.  A  further  observations  that  Can  be  made  on  the  results  of  the 
simulation  is  the  rapid  decay  of  the  resolved  turbulence  with  the  increase  in  the  width  of  the 
wake.  In  the  far  wake  only  larger  turbulence  structures  can  been  seen.  This  behavior  seems 
reasonable,  since  the  small  turbulence  structures  contain  significantly  less  energy,  hence  they  can 
only  last  for  shorter  time  compared  to  the  larger  eddies. 


Figure  4.2:  a)  Geometry  of  DTMD  5415,  b)  Grid  of  the  computational  domain:  162x50x66. 
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The  mean  kinetic  energy  was  calculated  from  the  predicted  velocity  field.  This  was 
accomplished  by  calculating  the  mean  velocity  on  several  cross-sections  in  the  flowfield  via 
summation  of  the  instantaneous  velocity  over  20000  timesteps.  Then,  these  mean  velocity  fields 
were  subtracted  from  the  instantaneous  velocity  field  on  those  cross-sections  to  calculate  the  root 
mean  square  velocity  fluctuations  and  hence  the  turbulent  kinetic  energy  field. 

Figure  4.8  shows  that  noticeable  changes  are  observed  in  the  streamwise  velocity  profiles 
in  both  axial  and  vertical  directions.  The  velocity  has  a  peak  below  the  free  surface,  which  gets 
weaker  in  the  axial  direction.  However,  the  location  of  peak  axial  velocity  defect  moves  deeper 
into  the  wake  indicating  a  plunging  effect  as  the  wake  develops  as  seen  in  Figures  4.10(a)- 
4.10(c).  The  turbulent  kinetic  energy  of  the  rms  velocity  fluctuations  in  Figure  4.9  shows  decay 
and  significant  redistribution.  Initially  the  turbulence  is  concentrated  more  near  the  free  surface. 
This  area  plunges  deeper  with  increasing  axial  distance  as  also  observed  in  Figures  4.11(a)- 
4.1 1(c).  This  phenomena  is  investigated  in  more  detail  by  Yavuz  et  al.,  2004.  Moreover,  it  can  be 
concluded  from  theses  figures  that  the  LES  predictions  clearly  depict  the  two  bilge  vortices, 
whereas  the  RANS  predictions  are  known  to  be  incapable  to  do  so. 

The  effects  of  different  numerical  schemes  and  SGS  models  were  also  investigated  in  this 
study.  The  resolved  turbulence  kinetic  energy  (TKE)  for  different  schemes  and  different  grid 
sizes  is  shown  in  Fig.  4.12.  The  RANS  solution  obtained  by  the  CFD-SHIP  code  from  IOWA  is 
also  shown  in  this  picture  for  comparison.  The  resolved  TKE  from  fine  grid  is  higher  than  that  of 
the  coarse  grid  as  expected.  Central  differencing  discretization  with  Smagorinsky  model  gives 
better  results  than  the  other  schemes  in  retrospect  to  the  total  TKE  calculated  from  RANS.  From 
this  figure  it  can  also  be  seen  that  there  is  no  significant  difference  in  the  results  when  the 
QUICK  scheme  is  used  with  or  without  SGS  model.  Moreover,  the  resolved  TKE  is  lower  than 
that  of  central  differencing  with  Smagorinsky  model.  This  means  that  the  grid  is  relatively  coarse 
and  the  QUICK  scheme  contains  even  higher  numerical  diffusion  than  the  Smagorinsky  model. 
For  more  detailed  information  about  the  comparisons  of  numerical  schemes  and  subgrid-scale 
models  the  reader  is  referred  to  Shi  et  al.  (2000b)  and  Shi  (2001).  One  uncertainty  in  the 
computations  is  the  presence  of  a  sinusoidal-like  distribution  of  TKE  in  the  near  wake.  It  may  be 
because  of  the  existing  surface  wave  (not  accounted  directly  here).  When  the  wave  descends 
towards  the  bottom  of  the  domain,  it  creates  a  constriction  with  flow  passing  through  a  small 
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area.  Thus  both  the  velocity  and  TKE  are  higher  at  this  region.  A  closer  look  at  the  wave  profile 
of  RANS  calculation  showed  that  the  peak  region  of  k-profile  did  indeed  correspond  to  the 
descending  wave  region.  It  is  interesting  that  our  calculations  also  provide  a  similar  trend 
although  no  surface  wave  profile  was  applied  in  our  LES.  This  indicates  that  some  wave 
information  may  be  present  implicitly  in  the  inflow  boundary  data. 

Figures  4.13  and  4.14  show  the  inlet  flow-field  as  computed  with  RANS  model  and 
produced  after  the  application  of  the  RFG  procedure,  respectively.  In  figures  4.15-4.19  the 
velocity  vectors  on  different  vertical  cross  sectional  planes  are  presented.  Clearly,  the  large  scale 
turbulent  eddies  (vortices)  are  captured.  As  can  be  seen  from  Fig.  4. 1 7  the  axial  development  of 
the  flow  shows  the  existence  of  unsteady  vortices  in  the  region  around  X/L  =  1 .  These  consist  of 
two  stable  large  bilge-type  vortices,  and  an  additional  smaller  vortices  farther  away  from  the 
axis.  The  axial  development  of  the  flow-field  causes  these  vortices  to  weaken  and  drift  in 
spanwise  direction  toward  the  outer  region  of  the  wake,  where  they  eventually  disappear  (Fig. 
4.19).  Although  there  are  no  bilge  vortexes  at  the  inlet  boundary,  these  structures  are  developed 
in  the  wake,  which  can  be  clearly  seen  in  Figs.  4.15-4.19.  This  result  gives  more  credit  to  our 
approach  of  inlet  boundary  specifications,  which  inherits  most  of  the  flow  characteristics  after 
the  ship  hull.  The  prediction  of  outward  flow  near  free  surface,  and  the  streamwise  evolution  of 
the  vortices,  without  largely  being  dissipated,  are  also  very  encouraging  indications  of  the 
success  of  the  present  LES  approach. 

The  vorticity  contours  on  different  vertical  planes  shown  in  Figs.  4.20  -  4.22  are 
indicative  of  the  degree  of  resolution  of  the  calculations.  These  structures  are  hard  to  see  in 
Figures  4.15-4.19  due  to  relatively  small  magnitude  of  the  velocity  vectors.  However,  to 
demonstrate  the  small  weak  turbulence  structures,  like  those  in  Fig.  4.18(a),  we  enlarged  the 
velocity  vectors  in  the  areas  containing  turbulence  structures  corresponding  to  those  in  the 
vorticity  contour  plot  and  depicted  them  in  Fig.  4.18(b).  A  corresponding  area  in  the  vorticity 
contour  plot  is  also  included  in  Fig.  4.22  for  reference  with  the  cross  section  taken  at  x/L  =  1 .2. 
Without  interpolation,  the  weak  vortices  can  not  be  seen  in  the  velocity  vector  plot  due  to  scale 
difference,  whereas  the  interpolated  contour  plots  do  show  small  structures. 

Figures  4.21  and  4.22  show  again  that  much  of  the  vorticity  is  concentrated  near  the  free 
surface,  and  there  are  two  large  counter  rotating  vortexes  on  two  sides  of  the  wake  in  addition  to 
several  smaller  ones  which  can  play  an  important  role  in  bubble  dynamics. 
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Figure  4.9:  Turbulent  kinetic  energy  variation  in  axial  and  vertical  directions. 
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Figure  4.10:  Axial  mean  velocity  contours  in  the  wake  of  the  ship  model  5512. 
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Figure  4.1 1:  Turbulent  kinetic  energy  contours  in  the  wake  of  the  ship  model  5512 
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Figure  4.12:  Comparison  of  resolved  turbulence  kinetic  energy  for  different  simulated  cases;  the 

wake  of  a  ship  model  5512 
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Figure  4.13:  RANS  solution  at  x/L  =  0.0 
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Figure  4.14:  RANS+RFG  solution  at  x/L  =  0.0 
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Figure  4.19:  Typical  instantaneous  velocity  vectors  on  y-t  plane  at  x/L  =  1 .4  in  the  wake  of  the 

ship  model  5512 


Figure  4.21:  Typical  instantaneous  vorticity  contours  (cox)  on  y-z  plane  at  x/L  =  0.6  in  the  wake  of 


the  ship  model  5512 
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Figure  4.22:  Typical  instantaneous  vorticity  contours  (cox)  on  y-z  plane  at  x/L  =  1 .2  in  the  wake  of 

the  ship  model  5512 
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4.4  Ship  Wake  Simulations  on  a  Circular  Track 

As  a  first  step  for  this  study,  a  standard  Smagorinsky  model  that  forms  the  basis  for  the 
advanced  models  has  been  applied  to  investigate  the  influence  of  the  Coriolis  and  centrifugal 
forces  on  turbulence  generation  on  the  ship  wake  flow  field. 

The  focus  of  this  study  is  the  wake  behind  the  Navy  DDG51  surface  ship,  which  is 
approximately  154  m  long  and  20  m  wide,  cruising  on  a  circular  track  (see  Figure  4.23).  The 
average  speed  of  the  ship  is  assumed  to  be  around  20  knots.  The  ship  model  data  used  for  this 
simulation  was  taken  from  the  data  of  DTMB  5415,  a  towing  tank  model  representing  a  modern 
naval  combatant,  DDG51  (Stern  et  ah,  2000).  As  stated  in  Shi  (2001)  and  Yavuz  et  al.  (2002), 
the  ship  hull  is  excluded  from  the  LES  calculations  due  to  computational  expenses  and  the 
computations  were  started  from  a  plane  aft  of  the  ship.  This  is  accomplished  using  the  RFG 
technique,  originally  developed  at  West  Virginia  University,  which  calls  for  a  time  averaged 
flow  field  at  the  inlet  data  plane.  Reynolds  Averaged  Navier  Stokes  (RANS)  calculations 
(Hyman,  2001)  (steady  state  RANS  calculations  around  the  ship  hull)  are  used  to  provide  the 
RFG  procedure  with  the  information  needed  on  the  inflow  boundary  (initial  data  plane)  located 
0.5L  after  the  body  in  the  wake.  In  other  words,  the  pseudo  random  flow  field  generated  by  the 
RFG  technique  is  added  to  the  mean  flow  of  the  RANS  simulations  in  order  to  provide  the 
boundary  condition  at  the  inlet  plane.  The  further  development  of  the  wake  flow  is  calculated  via 
LES  (non  steady  LES  of  the  wake).  LES  of  high  Reynolds  number  flows  with  complicated 
geometries  are  enabled  by  this  combined  approach  (Smirnov  et  al.,  2000).  Therefore,  the  effect 
of  the  ship  body  on  the  flow  field  is  embedded  in  the  mean  flow  prescribed  at  the  inflow  plane. 
However,  this  technique,  also  known  as  the  Initial  Data  Plane  (IDP)  approach  (Hyman,  1998, 
Paterson  et  al.,  1996)  can  introduce  considerable  errors  (Hyman,  1998),  as  a  matter  of  fact  that 
the  body  generating  the  wake  is  not  included  in  the  non-steady  LES  simulations,  It  should  be 
noted  that  it  is  possible  to  predict  turbulence  via  LES  technique  by  starting  from  a  mean  flow 
obtained  from  RANS.  However,  it  takes  a  very  long  time  for  a  turbulent  flow  to  develop  spatially 
and  temporally  without  any  initial  perturbation.  Also  for  some  engineering  applications,  it  may 
be  too  expensive  in  terms  of  computer  resources  and  programming  effort. 
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4.4.1  Conditions  and  Grids 

As  stated  above,  the  ship  model  DTMB  5415  is  a  5.72m  long  model  of  the  Navy  DDG51 
surface  ship  (Stem  et  al.,  2000).  A  computational  domain  of  1.75x0.15x1.0  (given  in  non 
dimensional  units  in  ship  lengths  and  it  starts  from  x/L=1.50,  where  x/L=l  is  the  end  of  the  ship 
model)  and  a  grid  of  66x33x66  (Coarse  Grid)  (Case-1)  in  x  (axial),  y  (vertical),  z  (transverse) 
directions  has  been  used  to  represent  the  near  wake  region.  This  grid  configuration  has  been 
selected  to  quickly  assess  the  effect  of  the  Coriolis  and  centrifugal  forces  on  the  turbulence 
characteristics.  It  should  be  noted  that  a  thorough  study  was  not  conducted  as  to  whether  the  side 
boundary  is  far  enough  away  from  the  wake.  In  the  IDP  plane,  the  RFG  method  is  used  in 
conjunction  with  the  RANS  calculations  (Hyman,  2001)  and  the  ship’s  stern  is  at  (0.5,  0,  -3.0)  in 
x-,  y-,  z-  directions  respectively.  The  core  region  of  the  numerical  grid  and  the  geometry  of  the 
ship  model  are  illustrated  in  Figure  4.24  (a)  &(c).  The  ship  turns  with  a  dimensionless  angular 
velocity  of  1/3  (Figure  4.24  (b))  with  a  radius  of  curvature  corresponding  to  three  dimensionless 
ship  lengths,  which  will  result  in  a  dimensionless  ship  velocity  of  1.  Only  a  30°  turn  was 
investigated.  Two  Reynolds  numbers  have  been  used  in  the  simulations.  One  is  based  on  the  real 
ship  length,  i.e.  1.5xl09  stated  as  Case-1  and  the  other  based  on  the  model  ship  length,  i.e. 
l.OxlO7  stated  as  Case-2.  Here,  Reynolds  number  similarity  is  assumed,  which  is  attained  by 
changing  the  laminar  viscosity  of  the  fluid  in  the  simulations. 

The  coordinate  system  used  is  with  respect  to  an  observer  on  the  ship.  The  Smagorinsky 
constant  used  in  the  eddy  viscosity  relation  is  calculated  to  be  0.042  from  scaling  with  the  non- 
turning  simulations  (Shi,  2001).  Moreover,  a  simulation  with  relatively  fine  grid  in  the  axial 
direction  has  been  conducted  to  investigate  the  grid  sensitivity  of  the  predictions.  This  grid 
consists  of  130x33x66  nodes  (Medium  Grid)  for  Case-1  in  x,  y,  and  z  directions,  respectively. 
Non-uniform  grids  were  used  in  this  study  in  both  x  and  z  directions  with  the  expansion  ratio  not 
exceeding  1.03.  The  length  scale  and  time  scale  used  in  RFG  are  calculated  from  x-k/s  and 
l=0.09k1,5/8  ,  non  dimensionalized  by  free  stream  velocity  and  ship  length.  A  comparison  of  the 
results  from  using  Smagorinsky  constants,  0.065  (selected  to  be  the  constant  for  the  non-turning 
ship  simulation)  and  0.042  has  been  presented  for  the  relatively  fine  simulations.  Then,  the  grid 
of  130x50x110  (Fine  Grid)  and  190x50x110  (Finest  Grid)  for  Case-2  have  been  used  to  analyze 
the  physics  of  the  wake  behind  a  turning  ship.  Finally,  to  investigate  the  effect  of  the  free 
surface,  another  study  has  been  conducted  with  fine  grids  for  the  model  ship  simulation.  The 
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time  step  is  0.001  for  all  grid  resolutions.  On  an  Intel  Pentium  4  3GHz  machine,  CPU-time  and 
memory  requirements  were  approximately  226  hours  for  otie  flow  through  time  (time  required 
for  the  flow  to  pass  in  the  calculation  domain,  x=1.75L)  ahd  491  MB  respectively,  for  the  fine 
grid  simulations  with  the  standard  Smagorinsky  model. 

On  the  other  hand,  Shi  (2001)  and  Shi  et  al.  (2006)  applied  LES  on  the  wake  flow  of  the 
ship  model  DTMB  5512  (Longo  et  al.,  1993,  Gui  et  al.,  1999,  Stem  and  Wilson,  2000).  A 
3.048m  long  unpropelled  model  of  a  modem  U.S.  Navy  fleet  ship,  Arleigh-Burke  class 
destroyer,  DDG51,  with  a  Reynolds  number  of  4.65x106  cruising  on  a  straight  track  has  been 
investigated.  The  computational  domain  was  1.5x0.3x0.6  (given  in  non  dimensional  units  in  ship 
length)  in  x- ,  y-,  z-directions,  respectively,  with  a  grid  size  of  162x50x66  and  322x50x66.  Non- 
uniform  grid  spacing  with  stretching  smaller  than  1 .03  was  used  in  both  x  and  y  directions. 


Figure  4.23.  Arleigh-Burke  class  destroyer  (DDG51) 


4.4.2  Boundary  Conditions 

For  all  simulations,  the  velocity  components  at  inflow  boundary  plane  are  set  equal  to 
those  calculated  via  RANS  simulations  by  Hyman  (2001)  (see  Figure  4.24.  (b))  plus  the 
fluemating  velocity  component  obtained  from  the  RFG  technique.  The  outflow  boundary  is 
assumed  to  be  a  free  gradient  boundary.  At  the  top  and  bottom  boundary  surfaces  a  symmetry 
boundary  condition  is  applied.  At  the  free  surface  (i.e.  the  tOp  boundary),  a  slip  is  allowed  in  x 
and  z  directions  but  the  velocity  component  normal  to  the  frde  surface  is  set  to  zero  (as  such  the 
free  surface  is  approximated  as  a  moving  flat  plane;  where  the  free  surface  is  completely  flat). 
The  boundaries  in  the  transverse  direction  are  treated  as  sfream  surfaces,  where  the  tangential 
velocities  are  specified  using  the  turning  ship  velocity  as  a  base  and  by  accounting  for  the  turning 
of  the  ship  via  a  rigid  body  motion,  such  as, 
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Qs  =  the  angular  velocity  of  the  system  rotation  due  to  the  turning  ship 
R*  =  the  radius  of  the  curvature,  R*  =  3.00L 

z*  =  the  distance  measured  from  the  centerline  in  the  positive  z  direction 


Figure  4.24.  Turning  ship  wake:  a)  The  geometry  b)  velocity  profile  specified  at  the  IDP  (top 
view)  c)  The  coordinate  space  system  and  numerical  grid  (Only  the  core  region  is  shown,  distances  are 

non-dimensionalized  with  ship  length,  L) 


4.4.3  Results  and  Discussions 

Figure  4.25  shows  the  instantaneous  vertical  velocity  contours  with  Coriolis  force  (Case- 
1)  at  different  y-z  planes  in  the  turning  ship  wake.  This  figure  has  been  included  to  help  the 
reader  to  visualize  the  cross  sectional  planes  taken  at  different  angles  (5°,  10°,  15°,  20°,  25°,  30°). 
Figure  4.26  shows  the  radial  velocity  profiles  at  two  different  y-z  planes  (5°,  10°)  for  fine  grid  in 
order  to  check  the  boundary  conditions  specified  in  the  transverse  direction.  The  figure  confirms 
that  the  LES  code  treats  the  boundaries  in  the  transverse  direction  as  stream  surfaces.  The 
unsteady  velocity  fluctuations  are  compared  for  cases  with  and  without  Coriolis  force  (Case-1)  in 
Figures  4.27  to  4.28.  The  figures  show  the  streamwise  velocity  components  versus 
nondimensional  time  at  an  angle  of  10  degrees  measured  from  the  inlet  plane  z.  From 
approximate  calculations,  the  resolved  frequency  of  the  velocity  fluctuations  is  obtained  to  be 
around  4.7  Hz  in  all  three  directions  with  the  Coriolis  force  and  3.5  Hz  without  the  Coriolis 
force.  This  indicates  that  the  Coriolis  force  tends  to  increase  turbulence  activity.  In  order  to 
compare  the  results  of  the  turning  and  non  turning  ship  wake  cases,  a  scaling  has  been  done  such 

that,  (C.A2)  =(CS A2)  .  From  Shi’s  (2001)  non  turning  ship  wake  calculations,  Cs  is 

nontum  turn 

set  to  be  0.065.  Using  the  above  scaling,  Cs  of  0.042  has  been  calculated  and  used  for  the 
medium  grid  simulations. 

For  fine  and  finest  grid  simulations,  similar  logic  has  been  carried  out  in  order  to  obtain 
Cs  values  in  the  Smagorinsky  model.  For  the  medium  grid,  the  resolved  frequency  is  around  9Hz 
for  Cs=0.042  and  7  Hz  for  Cs=0.065.  The  frequencies  are  approximately  calculated  from  the 
temporal  history  curves  of  all  three  directions,  such  that  one  over  the  difference  of  the  two  crest 
points  divided  by  the  time  difference  in  those  crest  points  gives  the  approximate  frequency.  If  the 
Smagorinsky  constant  is  smaller,  the  frequency  obtained  is  higher,  which  implies  that  more 
energetic  turbulent  fluctuations  are  captured.  This  is  an  indication  that  higher  frequencies  and 
smaller  turbulence  scales  are  captured  with  the  relatively  fine  grid  for  both  values  of  the 
Smagorinsky  constant.  Hence,  physical  intuition  and  simulations  give  similar  conclusions. 

Figures  4.30-4.32  present  the  velocity  vectors  at  vertical  cross-sections  at  an  angle  of  5° 
with  the  z-axis  Note  that  R  is  taken  to  be  minus  to  be  consistent  with  the  minus  z-  direction,  that 

is  calculated  from  R  =  yj(x  —  0.5)2  +  y2  .  The  coarse  grid  case  in  Figure  4.30  with  the  Coriolis 
force  shows  better-defined  turbulent  structures  or  eddies  than  the  case  without  the  Coriolis  force 
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in  Figure  4.31,  which  seems  to  be  more  diffusive.  Here,  oh#  could  argue  that  the  Coriolis  force 
helps  to  maintain  the  turbulent  structures  of  the  flow.  Moreover,  the  medium  grid  predictions 
with  the  Coriolis  force  presented  in  Figure  4.32  indicate  even  a  lesser  diffusion  (or  smearing)  for 
both  values  of  Smagorinsky  constant.  Here,  it  can  be  seen  that  the  vortical  structures  are  not 
penetrating  as  deeply  as  in  the  coarse  grid  case.  Experimental  observations  of  Matsubara  and 
Alffedsson  (1998)  have  shown  that  the  Coriolis  force  may  give  rise  to  instabilities  in  the  form  of 
longitudinal  vortices  which  supports  the  above  assertions. 


Figure  4.25.  Typical  instantaneous  vertical  velocity  contours  with  Coriolis  force  on  different  y-z 
plane  in  the  turning  ship  wake  (coarse  grid) 


Figure  4.26.  Radial  velocity  profiles  on  different  y-z  planes  in  the  turning  ship  wake  (fine  grid) 
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Figure  4.27.  Temporal  history  of  streamwise  velocity  components  at  10°with  the  z  axis,  and 


x=1.02,  y=-0.002  (coarse  grid),  a)  w  Coriolis  force  b)  w/o  Coriolis  force 


Figure  4.28.  Temporal  history  of  streamwise  velocity  components  at  10°with  the  z  axis,  and 


x=1.02,  y=-0,002  (medium  grid),  Cs=0.042 


Figure  4.29.  Temporal  history  of  streamwise  velocity  components  at  10°  with  the  z  axis,  and 


x=T.02,  y=-0.002  (medium  grid),  Cs=0.065 
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Figure  4.30.  The  velocity  vectors  on  a  vertical  plane  at  an  angle  of  5°  with  the  z-axis  with 

Coriolis  force  (coarse  grid) 


Figure  4.31.  The  velocity  vectors  on  a  vertical  plane  at  an  angle  of  5°  with  the  z-axis  without 


Coriolis  force  (coarse  grid) 


The  velocity  vectors  on  vertical  cross-sections  at  an  angle  of  25°  with  the  z-axis  are 
shown  in  Figures  4.33-4.34.  Although,  similar  structures  in  both  near  the  inlet  data  plane  and 
near  the  outlet  are  observed,  the  intensity  of  the  velocity  fluctuations  diminish  as  the  outflow 
plane  is  approached  due  to  grid  expansion.  Moreover,  it  is  seen  that  away  from  the  wake 
centerline,  the  strength  of  the  flow  structures  has  already  died  in  both  cases.  In  the  far 
downstream  portion  of  the  calculation  domain,  larger  turbulent  structures  are  seen  to  be  merging 
with  the  smaller  structures.  Still,  a  significant  level  of  vorticity  can  be  captured  in  both  cases. 
Overall,  it  is  seen  that  one  large  vortex,  presumably  the  one  originating  from  the  ship  hull 
persists  without  significant  dissipation. 

A  comparison  of  the  resolved  turbulence  kinetic  energy  is  shown  in  Figure  4.36.  In 
calculating  the  kinetic  energy  values,  the  samples  were  taken  to  be  1 0000  time  steps  from  2  flow 
through  times  data  (=35000  time  steps).  It  is  observed  that  due  to  the  coarse  grid  towards  the  end 
of  the  calculation  domain,  vorticity  (hence  turbulence)  decays  rapidly,  however  the  medium  grid 
predictions  have  a  very  gradual  decay  of  turbulent  kinetic  energy.  For  the  coarse  grid 
simulations,  turbulence  decays  more  rapidly  due  to  the  numerical  dissipation.  The  kinetic  energy 
values  obtained  from  the  case  with  Coriolis  force  is  slightly  higher  than  that  of  the  case  without 
Coriolis  force,  as  expected.  As  seen  from  the  figure,  the  adjustment  of  the  Smagorinsky  constant 
is  necessary,  as  the  kinetic  energy  values  obtained  with  Cs=0.042  are  almost  on  the  same  level  as 
the  ones  obtained  from  the  non-turning  ship  simulations(Shi,  2001).  These  observations  allude  to 
the  fact  that  the  Smagorinsky  constant  should  be  adjusted  for  streamline  curvature  effects.  The 
turbulence  intensity  specified  at  the  inlet  was  observed  to  decay  rapidly,  which  was  mainly  due 
to  the  nature  of  the  coarse  grid  used.  For  this  reason  a  medium  grid  simulation  has  been 
performed.  This  improved  the  turbulent  kinetic  energy  prediction  and  more  detailed  turbulence 
structures  were  captured. 

Similar  conclusions  could  also  be  drawn  from  the  vorticity  contours.  Calculations  showed 
that  the  highest  and  lowest  vorticity  magnitudes  are  more  pronounced  in  the  case  with  the 
Coriolis  force,  whereas  the  case  without  the  Coriolis  force  shows  a  more  diffused  or  smeared 
picture.  The  vorticity  results  again  indicate  that  the  Coriolis  force  keeps  the  vorticity 
concentrated  in  the  flow-field. 
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Figure  4.33.  The  velocity  vectors  on  a  vertical  plane  at  an  angle  of  25°  with  the  z-axis  with 

Coriolis  force  (coarse  grid) 


Figure  4.34.  The  velocity  vectors  on  a  vertical  plane  at  an  angle  of  25°  with  the  z-axis  without 

Coriolis  force  (coarse  grid) 


Figure  4.35.  The  velocity  vectors  on  a  vertical  plane  at  an  angle  of  25°  with  the  z-axis  (medium 

grid)  a)  Cs=0.042  b)  Cs=0.065 
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Figure  4.36.  The  comparison  of  the  resolved  turbulence  kinetic  energy  for  ship  cruising  on  a 

circular  track 
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Figure  4.38.  Velocity  vectors  at  x/L=0.65  a)  Non-turning  ship  (Shi  et  al.  2001)  (finest  grid)  b) 

Turning  ship  (finest  grid) 


Figure  4.39.  Velocity  vectors  at  x/L=0.65  for  turning  ship  wake  a)  Case-1;  Re=1.5xl09:  DDG51 
(finest  grid)  b)  Case-2;  Re=l. 0x107:  DTMB  5415  (finest  grid) 

The  goal  of  this  preliminary  study  was  to  assess  the  effects  of  the  Coriolis  and  centrifugal 
forces  on  the  vortical  structures  or  turbulence  characteristics  of  the  flow  in  the  wake  of  a  turning 
ship  using  the  large  eddy  simulation  technique.  The  eddies  resolved  by  LES  have  been  observed 
to  be  more  energetic  and  less  diffusive  when  the  Coriolis  force  was  included.  It  seems  as  if  this 
force  supplies  energy  to  the  large  turbulent  structures  and  thus  enhances  anisotropy.  The  vorticity 
contours  show  a  non-symmetric  wake  development  with  significant  stretching  in  the  radial 
direction  away  from  the  center  of  rotation.  This  is  also  seen  from  the  velocity  vectors  comparison 
of  non-turning  and  turning  ship  studies  in  Figures  4.37  &  4.38.  As  these  figures  show  the  non- 
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turning  ship  case  has  a  symmetric  wake  with  respect  to  the  wake  centerline,  whereas  in  the 
turning  ship  case,  there  is  a  significant  flow  present  in  the  radial  direction  that  is  believed  to  be 
caused  in  part  by  the  centrifugal  force,  arising  from  turning  of  the  ship.  Here,  the  simulation  of 
the  model  ship  with  grids  of  130x50x110  and  190x50x110  has  been  studied  to  compare  the 
results  with  the  non-turning  ship  wake. 

A  study  of  the  effect  of  Reynolds  number  for  Case-2  and  Case-1  using  190x50x1 10  grids 
is  shown  in  Figure  4.39.  This  comparison  can  be  done  since  all  the  values  are 
nondimensionalized  with  respect  to  their  ship  velocities,  either  real  ship  or  the  ship  model.  The 
velocity  vectors  are  at  the  same  location,  x/L=0.65  and  standard  Smagorinsky  model  has  been 
used.  When  the  Reynolds  number  is  low,  the  turbulence  structures  are  observed  to  be  more 
visible,  well  defined  and  more  diffusive.  However  for  high  Reynolds  number  flows,  there  are 
certain  unorganized  structures  seen  in  the  flow  field  and  the  flow  structures  seem  to  be  less 
diffused.  As  the  wake  develops,  the  flow  structures  disappear  rapidly  for  Case-1;  however,  for 
Case-2,  the  flow  structures  are  observed  to  be  still  well  defined.  The  unsteady  velocity 
fluctuations  are  also  compared  for  Case-1  and  Case-2  in  Figures  5.40.  The  resolved  frequency  of 
the  velocity  fluctuations  is  estimated  to  be  around  40  Hz  in  all  three  directions  for  Case-1  and  30 
Hz  for  Case-2.  Here,  30  Hz  and  40  Hz  seem  to  correspond  to  large  bilge  vortex  passage.  If  the 
Re  is  higher,  the  frequency  obtained  is  higher,  which  is  an  indication  that  smaller  turbulent 
structures  are  captured  with  high  Reynolds  number.  The  period  (hence  frequency)  can  also  be 
approximated  from  the  roughly  eddy  turnover,  which  is  T=largest  eddy  size/ship  velocity  where 
at  10°,  the  largest  eddy  size  is  approximately  0.04  from  Figure  4.45  (b)  and  the  ship  velocity  is  1, 
therefore  T=0. 04/1  =0.04,  f=l/T=25Hz,  which  is  close  to  30Hz. 

Although  the  classical  Smagorinsky  SGS  model  is  not  that  suitable  for  complex  flows  as 
it  uses  a  constant  eddy  viscosity  coefficient  for  the  entire  domain,  this  study  has  shown  that  it  can 
be  used  as  a  SGS  to  predict  the  flow  dynamics  of  the  wake  behind  a  turning  ship  reasonably 
well. 

4.4.4  Properties  of  Turbulent  Ship  Wakes 

The  mean  inflow  boundary  data  was  obtained  by  slicing  the  RANS  solution  at  x/L=T.5 
plane  for  turning  ship  wake  and  x/L=l  .05  for  the  non  turning  ship  wake  and  then  interpolated  to 
the  inlet  plane  of  the  computational  domain.  The  mean  axial  velocities  for  both  non  turning  and 
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turning  ship  wake  after  the  interpolation  are  presented  in  Figure  4.41  (a)  and  (b).  To  see  the 
effect  of  the  LES  wake  calculation  only,  the  axial  velocity  contours  are  adjusted  due  to  rotation 
and  the  mean  axial  velocity  adjusted  for  rotation  at  the  IDP  (x/L=1.5)  for  turning  ship  wake  is 
given  in  Figure  4.41  (c).  By  using  RFG  and  the  mean  flows  on  these  planes,  the  turbulent  inflow 
boundary  for  both  ship  simulations  was  reconstructed.  Figure  4.42  presents  the  pseudo  random 
flow  field  generated  by  the  RFG  on  the  IDP  for  turning  ship  wake.  From  the  axial  velocity 
contours  in  Figure  4.43(a),  two  stable  large  bilge  type  vortices  are  observed  for  the  non-turning 
ship  wake  simulation  using  the  standard  Smagorinsky  model  (Re=4.65xl06)  at  x/L=0.65.  Two 
small  side  vortex  pairs  are  observed  away  from  the  center  of  the  wake.  These  vortex  pairs  drift  in 
the  spanwise  direction,  get  weaker  and  eventually  disappear.  From  the  axial  velocity  contours 
adjusted  for  rotation  at  x7L=0.65  for  the  turning  ship  wake  in  Figure  4.43(b),  one  large  bilge 
vortex,  that  moves  downwards,  is  observed  and  is  probably  the  result  of  merger  of  the  two 
vortices  (bilge  vortices  of  opposite  rotation)  under  the  action  of  Coriolis  and  centrifugal  forces 
for  the  turning  ship  wake  simulation.  There  is  a  smaller  circulation  region  (a  side  vortex)  on  the 
outer  rim  of  the  wake.  The  streamwise  flow  field  causes  the  side  vortex  to  weaken  and  the  wake 
decays  in  the  outer  region  of  the  near  wake  similar  to  non-turning  ship.  The  Coriolis  force  seems 
to  generate  more  energetic  and  less  diffusive  eddies,  hence  it  seems  to  increase  kinetic  energy 
content  of  the  wake  (Yavuz  et  al.,  2002). 


Figure  4.40.  Temporal  history  of  streamwise  velocity  components  at  10°  (x=1.02,  y=-0.001  and  z=-3.18) 

(finest  grid)  and  Cs=0.052  a)  Case-1  b)  Case-2 
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Figure  4.41  a)  The  mean  axial  velocity  contour  at  the  inlet  data  plane  (IDP)  for  non-turning  ship 
(Shi  et  al.  2001)  (finest  grid)  b)  The  mean  axial  velocity  contours  at  the  IDP  for  turning  ship  wake  (finest 
grid)  c)  The  mean  axial  velocity  contours  after  subtracting  solid  body  rotation  contribution  at  the  IDP  for 

turning  ship  wake  (finest  grid) 
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Figure  4.42.  a)  The  axial  flow  field  provided  by  RFG  at  the  inlet  plane  (IDP)  for  turning  ship 
wake  b)  Enlarged  view  of  A  (finest  grid) 


Figure  4.43.  a)  Axial  velocity  contour  at  x/L=0.65  for  non-tilfhing  ship  (Shi  et  al.  2001)  (finest  grid) 
b)  Axial  velocity  contours  after  subtracting  solid  body  rotation  contribution  at  x’/L=0.65  for  turning  ship 

wake  (finest  grid) 


\ 
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Figure  4.44  shows  a  comparison  of  the  predicted  velocity  contours  with  the  macro  wake 
measurements  at  x/L=1.20  by  Hoekstra  &  Ligtelijn  (1991)  for  the  ship  model  No.  5452 
(Re~lxl07)  with  the  non-turning  ship  wake  at  the  same  location.  In  both  studies,  minimum  axial 
velocity  occurs  near  the  free  surface  of  the  center  of  the  wake,  as  seen  in  Figure  4.44.  Both  axial 
velocity  contours  resemble  an  upside  down  rimmed  hat  form.  Two  side  lobes  close  to  the  free 
surface  and  a  central  lobe  are  observed.  As  stated  by  Hoekstra  &  Ligtelijn  (1991),  these  contour 
forms  are  very  much  alike  for  all  ship  hulls.  Figure  4.44  indicates  that  from  the  straight  ship 
wake  simulation,  similar  physics  are  obtained  when  compared  with  the  macro  wake 
measurements  at  the  same  location.  Moreover,  the  extent  of  axial  turbulence  intensities  is  also  in 
reasonably  good  agreement  with  measurements,  as  shown  in  Table  4.1.  Figure  4.45  shows  the 
axial  development  of  the  wake  of  a  turning  ship.  The  bulk  movement  of  the  large  bilge  vortex 
can  be  seen  clearly  through  the  wake  as  it  moves  under  the  action  of  the  Coriolis  and  centrifugal 
forces.  The  comparison  of  the  energy  spectra  of  the  velocity  fluctuations  for  the  two  wakes  is 
shown  in  Figure  4.46.  It  is  seen  that  the  turning  ship  wake  has  more  energetic  fluctuating  eddies 
(the  values  are  almost  10  times  higher)  as  compared  to  the  non-turning  ship  at  the  same  location. 
This  finding  is  expected  as  the  initial  value  of  the  kinetic  energy  for  the  non-turning  ship  wake  is 
almost  8  times  smaller  than  the  turning  ship  wake.  For  the  turning  ship  wake,  the  smallest  eddy 
turnover  has  been  calculated  to  see  whether  it  is  correlated  or  not  with  the  time  step  used  in  the 
simulation;  At=lxl0'4.  The  smallest  eddy  turnover  (numerical  or  discretization  turbulence)  can 
be  defined  as,  the  eddy  size  over  the  ship  velocity,  where  the  resolved  eddy  size  is 

2^AXjAx2Ax3  and  the  ship  velocity  is  1.  At  10°,  Axi=4. 523x1  O'3,  Ax2=2.779xl0'3, 

Ax3=4.07xl0'3,  which  yields  an  eddy  turnover  time  of  4.52x1  O'4,  and  a  frequency,  f,  equal  to 
2520  Hz.  Since  this  frequency  is  much  larger  than  the  estimated  frequency  of  30  Hz,  there  seems 
to  be  no  correlation  between  them.  It  should  also  be  noted  that  the  shape  of  spectra  for  the 
turning  ship  case  is  quite  different  from  the  straight  track  case  (Fig.  4.46).  In  addition  to  these, 
the  wake  spreading  or  wake  width  for  the  non-turning  ship,  which  is  obtained  to  be  roughly 
w~x1/4  is  consistent  with  Buller  &  Tunaley  (1989)’s  measurements.  Milgram  et  al.  (1993)  and 
Hoekstra  &  Ligtelijn  (1991)  found  w~xl/5.  The  spreading  rate  of  the  turning  ship  wake  can’t  be 
easily  observed  like  in  the  non-turning  ship.  This  is  also  seen  from  the  predicted  vertical  vorticity 
contours  in  Figure  4.25.  It  may  be  due  to  the  grid  coarsening  towards  the  end  of  the  calculation 
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domain.  Moreover,  the  “inboard”  side  of  wake  has  a  much  sharper  edge  while  the  “outboard” 
side  is  more  diffusive.  This  might  be  due  to  instability  in  the  “outboard”  side. 


for  any  ship  hull 

-  maximum  value  of  root  mean  square  fluctuations 

location 

experimental 

numerical 

(m/s) 

(values  obtained  *US)  (m/s) 

x/L=0.25 

0.106±  0.02 

0.083 

x/L=0.6 

0.067  ±0.01 

0.069 

x/L=1.0 

0.047  ±0.01 

0.057 

Table  4.1 .  Comparison  of  the  maximum  values  of  the  root  mean  square  fluctuations  for  some 
locations  in  units  of  (m/s)  for  non-turning  ship 


Figure  4.44.  a)  Predicted  velocity  contours  at  x/L=l  .20  for  non-ttifhing  ship-(Shi  et  al.  2001)  (finest  grid) 
b)  Macro  wake  measurements  at  x/L=l  .20  by  Hoekstra  &  Ligtelijn  (1991 )  for  ship  model  No.  5452 


(Re~lxl07) 
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Figure  4.45.  Velocity  vectors  for  turning  ship  using  standard  Smagorinsky  model  (finest  grid)  on 
a  vertical  plane  at  an  angle  of  a)  5°  (x’/L=0.8)  b)  10°  (x’/L=1.02)  c)  20°  (x7L=1.55)  d)  30°  (x7L=1.75) 
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Figure  4.46.  Energy  spectra  of  the  velocity  fluctuatioiis  a)  Non-turning  ship  (Shi  et  al.  2001) 
(finest  grid)  b)  Turning  ship  (finest  grid) 


Figure  4.47.  Predicted  vertical  vorticity  (o>y)  for  a)  Non-turning  ship  (Shi  et  al.  2001)  (finest  grid),  b) 

Turning  ship  wake  (finest  grid)  (Case  2) 
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5  BUBBLE  DYNAMICS 


Prediction  of  bubble  distribution  in  ship  wakes  is  an  important  area  of  naval  research 
primarily  related  to  ship  detectability  from  air-bom  platforms.  Ship-wakes  can  be  located  by 
optical  means  using  detection  of  light  scattered  from  air  bubbles  entrained  in  the  wake.  On  the 
other  hand,  experimental  data  on  bubble  distributions  are  usually  obtained  by  acoustical 
methods.  Thus,  a  complete  analysis  of  this  problem  requires  complex  multi-physics  modeling 
that  would  account  for  different  physical  processes.  The  perspective  computational  tools  will 
have  to  be  based  on  multi-phase  fluid-dynamics  solvers,  capable  of  accurate  turbulence 
predictions,  and  with  potential  for  building  composite  models. 

Considering  relatively  small  volume  fractions  of  bubble  phase  in  the  ship  wakes,  the 
method  of  Lagrangian  particle  dynamics  (LPD)  is  most  appropriate  for  computing  bubble 
distributions.  It  has  advantages  over  the  Eulerian  two-fluid  methods  in  the  cases  of  dilute 
suspensions  or  when  large  concentration  variability  of  the  discrete  phase  is  present  [Elghobashi; 
Crowe].  This  situation  is  true  for  bubbly  wake  flows,  such  as  those  in  ship-wakes,  where  bubbles 
may  experience  preferential  concentration  and  clustering  effects  in  the  near  wake  region  but  are 
rather  dilute  in  the  far  wake.  The  LPD  method  is  also  conceptually  simple  and  methodologically 
robust  providing  a  numerically  stable  statistical  model  for  evaluating  dispersed  phase 
concentration  and  particle/bubble  distribution  functions.  The  algorithm  for  particle  tracking  and 
population  dynamics  developed  earlier  by  the  authors  demonstrated  the  ability  to  efficiently 
simulate  large  populations  of  particles  including  coalescence  effects  with  even  modest  computer 
resources  [Shi  et  al.  2000a&b;  Smirnov  and  Celik,  2000;  Smirnov  et  al.,  2000]. 

Simulations  of  bubbles  in  turbulent  shear  layers  require  accurate  representation  of  the 
fluctuating  flow-field  that  governs  bubble  dynamics.  When  conventional  RANS  (Reynolds 
Averaged  Navier-Stokes)  models  are  used  this  accuracy  is  often  lost  or  comes  at  a  cost  of 
empirical  physical  sub-models  for  turbulence.  On  the  other  hand,  direct  numerical  simulations 
(DNS)  or  LES  require  no,  or  a  relatively  simple  subgrid-scale  (SGS)  model.  LES  has  the 
advantage  over  DNS  of  being  able  to  handle  flows  with  higher  Reynolds  numbers  and  still 
accurately  reproduce  large-scale  dynamics  at  a  much  reduced  cost. 
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The  majority  of  simulations  of  turbulent  bubbly  shear  layers  using  Eulerian-  Lagrangian 
approach  is  done  using  RANS  methods  for  high  Reynolds  number  flows  [Joia,  1997;  Murai, 
2000]  or  DNS-type  simulations  for  relatively  low  Reynolds  number  flows  [Elghobashi  and 
Lasheras,  1996;  Ruetsch  and  Meiburg,  1994;  Okawa  et  al.,  2001;  Murai  et  al.,  2001],  Many 
commercial  CFD  codes  now  have  LES  algorithms  with  Lagrangian  tracking  techniques  for 
turbulent  flows.  However,  for  the  purpose  of  a  complex  multi-physics  analysis  of  turbulent 
wakes,  including  the  multi-phase  and  acoustical  effects,  the  authors  considered  it  important  to 
use  a  research-grade  LES  code.  This  will  enable  us  to  build  composite  solvers  including  fluid, 
acoustic  and  bubble  interactions.  Thus,  the  LES-based  approach  pursued  in  this  study  provides 
the  basis  for  multi-physics  LES  of  large-Reynolds  number  flows. 

One  of  perceived  disadvantages  of  the  Lagrangian  approach  is  the  requirement  to 
simulate  a  relatively  large  number  of  particle  trajectories  to  obtain  the  needed  statistics.  This 
study  elaborates  that  this  perception  is  exaggerated  as  long  as  the  flow  can  be  categorized  as 
nominally  dilute.  Particularly,  we  made  simple  estimates  (Sec.  5.4.2)  to  determine  under  what 
conditions  the  assumption  of  dilute  two  phase  flow  would  be  valid  for  calculations  of  turbulent 
bubbly  flow,  such  as  those  in  ship-wakes. 

5.1  Formulation 

The  LPD  algorithm  for  particle  tracking  and  population  dynamics  developed  by  the 
authors  demonstrated  the  ability  to  efficiently  simulate  large  populations  of  particles  including 
coalescence  effects  with  even  modest  computer  resources  [Shi  et  al.  2000;  Smirnov  and  Celik, 
2000;  Smirnov  et  al.,  2000;  Shi,  2000].  In  this  study  the  LPD  algorithm  was  combined  with  the 
Large-Eddy  simulation  approach  [Smirnov  et  al.,  2001;  Piomelli,  1999]  to  compute  the 
distribution  of  bubbly  phase  in  the  near-wake  flow  of  a  ship-model.  The  LES  technique 
[Piomelli,  1999]  enhanced  with  the  RFG  procedure  [Smirnov  et  al.,  2001;  Smirnov  et  al.,  2000; 
Shi,  et  al.,  2000]  was  implemented  into  the  LES  method  to  enable  the  appropriate  representation 
of  initial/inlet  turbulence  conditions  and,  as  a  subgrid  scale  for  particle  dynamics  solver.  It 
should  be  noted  that  the  RFG  approach  is  essentially  a  simplified  spectral  method,  and  as  such  it 
has  a  potential  to  be  used  as  a  part  of  an  acoustic  solver  in  the  appropriate  sub-model.  This  may 
be  important  for  validation  of  the  simulation  results  on  the  data  produced  in  acoustical 
experiments  [Trevorrow  et  al.,  1994], 
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The  combined  LES/RFG  method  was  validated  on  the  cases  of  turbulent  mixing  layer, 
and  used  later  to  simulate  the  wake  of  a  Navy  ship  model  541 54  [Smirnov  et  al.,  2001].  The 
fluctuating  velocities  at  the  inflow  boundary  provided  by  RFG  were  generated  from  the  data 
obtained  in  prior  RANS  calculations  based  on  k-e  turbulence  model  [Larreteguy,  1999], 

To  compute  bubble  motion  the  following  equations  were  adapted  from  Sridhar  and  Katz 
[Sridhar  and  Katz,  1995] 

~r~ =  A0  +  A,  +  Ad  +  A,  (5.1) 

at 

where  the  A  vectors  on  the  right-hand  side  represent  accelerations  due  to  added  mass 
( Aa ),  buoyancy  (  Aa  ),  drag  (Ad)  and  lift  ( A, )  given  by  the  following  expressions 
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(5.2) 


Relations  (5.2)  were  obtained  for  spherical  bubbles  and  by  considering  the  limit  of  a 
small  air-to-water  density  ratio  pjpv  «0.  The  coefficients  of  added  mass  ( Ca )  and  buoyancy, 
( Cb )  were  found  to  be  Ca  =3.0  and  Cb  =  2.0.  The  coefficients  of  drag  (Cd)  and  lift  (C;)  are 
themselves  empirical  functions  of  Ure/  and  bubble  diameter,  which  were  obtained  from 
experimental  measurements  of  bubbles  dynamics  in  turbulent  vortexes.  The  Bassett  term 
involving  the  history  integral  was  neglected  in  this  case,  following  the  conclusions  of  Sridhar 
and  Katz  [Sridhar  and  Katz,  1995]  on  small  contribution  of  this  term  compared  to  the  buoyancy 
term  for  a  typical  wake  flow  turbulence.  The  typical  bubble  sizes  used  in  the  experiments  were 
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500-800/i  m ,  and  the  Reynolds  numbers,  based  on  relative  velocity  between  20  and  80 ,  and 
vorticity  between  2  and  20U1 . 

The  set  of  equations  (5.1),  (5.2)  was  selected,  because  they  were  obtained  from  the 
experimental  data  on  bubble  dynamics  in  vortex  flows,  similar  to  those  of  turbulent  bubbly  wake 
flows.  Since  the  parameters  of  a  single  bubble  dynamics  in  simple  vortexes  were  directly 
measured,  these  data  represent  the  features  of  pure  bubble  dynamics  not  affected  by  the  random 
effects  of  turbulence,  as  may  be  the  case  with  some  experimental  data  collected  from  stochastic 
bubble  ensembles.  It  should  be  noted  that  there  are  other  alternatives  to  the  approach  selected 
here.  For  example,  the  classical  analytical  formulation  of  Auton,  [Auton,  1987;  Auton  et  al., 
1988]  and  it’s  extension  to  the  viscous  flow  [Auton  et  al.,  1988]  is  likewise  "non-polluted"  by 
stochastic  effects  of  turbulence,  and  this  form  of  the  lift  force  has  been  successfully  used  to 
predict  a  variety  of  bubbly  flows.  Generally,  the  question  of  which  expression  is  optimal  for  a 
given  bubbly  flow  is  still  an  open  research  issue.  The  authors  used  the  results  of  Sridhar  and 
Katz,  because  they  appeared  to  be  most  appropriate  for  representing  bubble  dynamics  in  a  vortex 
and  a  suitable  choice  for  LES.  However,  it  is  possible  that  similar  results  could  be  obtained  using 
another  set  of  equations. 

Equation  (5.1)  was  discretized  using  a  second-order  Runge-Kutta  time  stepping  scheme. 
In  the  mixing-layer  case  (Sec. 5. 3)  bubble  tracking  throughout  the  computational  domain  was 
simplified  by  using  a  uniform  Cartesian  grid,  so  that  the  index  of  the  cell  containing  the  bubble 
could  easily  be  obtained  using  a  simple  division  by  modulus  operation.  In  the  case  of  a  ship  wake 
(Sec. 5. 4)  this  restriction  was  lifted  and  a  more  general  grid-independent  tracking  algorithm 
(Sec. 5. 2)  was  used.  Fluid  velocities  were  interpolated  to  the  location  of  a  bubble  using  tri-linear 
interpolation  formula. 

In  a  joint  LES/LPD/RFG  approach  the  flow  solver  and  the  particle  solver  can  use 
different  time-stepping  schemes  with  independent  selection  of  time-steps.  This  was  necessary 
due  to  the  difference  in  turbulence  time-scales  and  bubbles  relaxation  times.  Usually  sub-cycling 
of  bubble  iterations  inside  the  flow  iterations  is  required  to  reach  a  stable  solution. 

To  estimate  the  dissolution  of  bubbles  of  different  diameters  at  different  depths  we  used 
the  bubble  dissolution  equation  as  formulated  in  [Rightley  and  Lasheras,  2000]  and  corrected  in 
[Hyman  1994;  Carrica  et  al.,  1998]. 
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where  H  is  Henry’s  constant,  PaIm  is  the  atmospheric  pressure,  g  is  the  gravity 

acceleration  constant,  and  a  is  a  surface  tension  coefficient. 

The  current  model  does  not  include  bubble  breakup.  This  can  be  justified  for  the  mixing 
layer  and  ship-wake  cases,  considered  below,  on  the  grounds  of  rather  small  turbulence  energies 
and  bubble  sizes.  Indeed,  if  we  use  the  conditional  breakup  probability  for  a  bubble  hit  by  an 
eddy  [Luo  and  Svendsen  1996  (Eqn.  26)]: 

12  cf<y  s 

",2/3  j5/3  ell/3 

PPcE  d  4  ) 

We  can  get  an  estimate  for  largest  bubbles  of  sizes  around  1  mm5.  The  breakup 
probability  in  this  case  is  around  1 0"5 ,  and  falls  rapidly  for  smaller  bubble  diameters,  d ,  as 
<sT5/3 ,  which  justifies  neglecting  their  breakup  effects. 


5.2  Implementation 

There  were  two  tasks  to  be  solved  for  an  efficient  implementation  of  an  LPD  procedure: 
(1)  particle  tracking  and  (2)  particle  population  dynamics.  The  first  task  is  related  to  the  problem 
of  identifying  the  grid-cell  where  the  particle  is  located,  and  interpolating  the  relevant  grid-based 
variables  to  the  particle  location.  The  second  task  is  related  to  the  problem  of  creating  and 
destroying  particles  as  the  need  arises  due  to  the  processes  of  particle  injection,  disappearance 
through  the  free  surface  and  crossing  the  boundaries  of  the  computational  domain.  Since  both 

5We  used:  cf  =  0.26 ,  a  -  0.0741  kg/s1 ,  (3  -  2.0,  pc  =1000  kg/nd ,  e  =1.0  m2/s 3 ,  £  =  1 
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particle  tracking  and  population  control  are  frequently  occurring  operations,  it  would  be  highly 
desirable  to  solve  each  task  in  a  manner  that  would  avoid  expensive  loops  through  the  particles 
or  grid-node  arrays. 

Particle  tracking 

When  the  particle  tracking  is  done  on  structured  grids,  the  location  of  the  grid-cell 
containing  the  particle  is  a  straightforward  task  reduced  to  a  division-by-modulus,  multiplication 
and  addition  operations.  When  non-uniform  or  unstructured  grids  are  used  the  procedure  is  not  as 
simple.  In  this  case  there  should  be  a  way  to  identify  if  the  current  location  of  the  particle  is 
inside  the  same  cell  it  was  inhabiting  in  the  previous  iteration  or  whether  the  particle  has  moved 
into  a  neighboring  cell. 

The  present  LPD  procedure  accomplished  this  task  by  solving  the  inclusion  problem  in  a 
polyhedron.  This  is  a  standard  algorithm  in  the  area  of  computational  geometry  [O'Rourke, 
1998],  where  the  given  point  is  tested  against  each  face  of  a  polyhedron  and  it  is  decided  if  the 
point  is  on  the  inside  or  the  outside  of  the  face  relative  to  the  center  of  the  polyhedron.  In  the  end 
the  point  is  either  found  to  be  inside  the  polyhedron  or  the  search  is  repeated  in  the  neighboring 
polyhedron  lying  across  the  first  face,  which  failed  the  test.  Identification  of  the  neighboring 
cells  is  a  simple  matter  on  structured  grids,  whereas  on  unstructured  ones  cell-neighbor 
connectivity  information  should  be  used6. 

Population  Dynamics. 

For  an  efficient  implementation  of  particle  population  control  algorithm,  it  would  be 
desirable  to  avoid  both  looping  trough  the  particles  and  dynamic  memory 
allocations/deallocations  in  cases  when  new  particles  are  created  or  destroyed.  It  can  be  achieved 
by  using  linked  lists.  In  the  present  scheme  an  array  of  the  size  equal  to  the  maximum  number  of 
particles  is  initially  allocated.  Then  the  pointer  system  is  setup  to  link  all  the  particles  together  in 
a  ring  fashion,  separating  them  into  active  and  dead  sublists  (Fig.5.1(a)).  Whenever  a  particle  is 


6In  parallel  implementation  the  neighbor  connectivity  should  cross  domain  boundaries 
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created  or  destroyed  it  takes  only  six  pointer  assignment  operations  to  move  it  between  the  active 
and  dead  sublists.  Figure  5.1(b)  shows  an  example  of  destruction  of  particle  B  by  the  appropriate 
reassignment  of  pointers  between  the  old  neighbor-particles  A,C  and  the  new  ones:  D,E. 

Even  though  the  limitation  on  the  maximum  number  of  particles  in  the  domain  may  seem 
restrictive,  it  is  a  compromise,  allowing  to  avoid  dynamic  allocation/  deallocation  of  memory 
space  for  each  newly  created  or  destroyed  particle.  Since  the  latter  operations  can  be  rather 
frequent  this  procedure  seems  reasonable  from  efficiency  standpoint. 


(a)  Particle  ring  (b)  Particle  destruction 

Figure  5.1  Dynamic  memory  allocation  of  particle  storage 


5.3  Validation  on  a  Turbulent  Mixing  Layer 

The  validation  of  the  LES/LPD  approach  was  done  on  experimental  data  of  a  mixing 
layer  [Rightley  and  Lasheras,  2000]  (Fig.  5.2).  A  detailed  description  of  the  validation  study 
including  the  results  for  the  flow-field  distribution  is  presented  in  [Smirnov  et  al.,  2001],  Elere 
we  present  additional  data  on  the  influence  of  the  effect  that  the  RFG  sub-grid  model  had  on  the 
results. 

In  the  experiments  the  mixing  layer  was  generated  by  two  separate  parallel  flows  with 
different  incoming  velocities.  A  thin  flat  plate  of  0.1 5m  length,  0.003m  height,  and  0.2m  width 
(the  whole  span-wise  extent),  was  mounted  in  the  middle  of  the  inlet  plane.  The  average  velocity 
of  the  lower  half  flow  was  0.28 m/s ,  while  the  upper  half  value  was  0.07 m/s .  Bubbles  of  40  /um 
in  diameter  were  carried  in  from  the  lower  half  of  the  channel.  For  further  details  of  the 
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experimental  installation  and  computational  setup  we  refer  to  the  relevant  sources  [Rightley  and 
Lasheras,  2000,15;  Smirnov  et  al.,  2001]. 


Figure  5.2  Mixing  layer  experimental  setup 


In  bubble  simulations  two  runs  were  conducted:  with  and  without  the  RFG  subgrid-scale 
model.  In  both  cases  a  column  of  100  bubbles  was  injected  at  the  inlet  of  the  domain  with  time 
interval  of  27  milliseconds  for  the  duration  of  12.7  seconds,  resulting  in  47,000  bubbles  that 
entered  the  domain  during  the  course  of  the  run.  Figure  5.3  shows  a  typical  bubble  distribution 
computed  by  LES/LPD  algorithm  compared  with  those  observed  in  experiments.  The  bubble 
cloud  can  be  seen  to  be  entrained  by  the  fluid  entering  the  mixing  region  from  the  high  speed 
side  into  the  cores  of  the  coherent  vortex  structures  present  in  the  mixing  region.  The  pictures 
represent  an  instantaneous  distribution,  which  is  different  for  any  given  time,  but  shows  common 
statistics  and  similar  dynamical  features. 

Spatial  histograms  of  bubble  distributions  were  obtained  by  counting  all  bubbles  passing 
the  cells  of  the  3D-histogram  and  accumulating  the  statistics.  Normalizing  the  histogram  data  by 
the  total  number  of  bubbles  injected,  gave  bubble  probability  density  functions  and 
concentrations.  The  statistical  error,  cr ,  in  the  number  of  bubble  counts,  n ,  for  each  slot  of  the 
histogram  can  be  estimated  from  binomial  distribution  as  cr=[n(l-o/JV)]1/2 ,  where  N  is  the 
maximum  number  of  bubbles  injected  during  the  simulation.  This  will  add  to  the  uncertainty  in 
the  layer  thickness  calculations  for  larger  X ,  and  smaller  bubble  concentrations. 

Figure  5.4  shows  the  surface  of  the  histogram  and  the  contour  plot  of  the  constant 
concentration  levels.  Figure  5.5  shows  the  growth  of  the  mixing  layer  thickness  in  the  stream- 
wise  direction  (X)  with  and  without  the  RFG  model.  Over-prediction  of  layer  thickness  at  small 
X  is  due  to  the  finite  resolution  of  the  histogram.  The  inclusion  of  the  RFG  model  improves  the 
predictions  at  greater  X  . 
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Figure  5.5  also  shows  predicted  and  measured  development  of  the  mixing  layer  thickness. 
The  overall  agreement  is  good  except  at  large  axial  distances  where  a  smaller  rate  of  growth  of 
the  mixing  layer  in  the  stream-wise  direction  is  observed  when  no  subgrid  scale  model  was 
applied.  Although  the  inclusion  of  the  RFG  model  [Smirnov  et  al.,  2001]  improved  the 
predictions  of  the  boundary  layer  growth  at  higher  axial  asymptotics,  it  also  increased  the 
uncertainty  in  the  boundary  layer  thickness  because  of  the  higher  dispersion  in  bubble 
distribution. 

Figure  5.6  provides  comparison  of  predicted  bubble  concentrations  against  experimental 
data.  The  agreement  is  reasonably  good  especially  at  the  beginning  of  the  boundary  layer.  Some 
irregularity  of  computed  data  further  downstream  is  most  probably  due  to  statistical  uncertainty 
of  the  sample  as  bubble  concentration  becomes  more  dilute. 

It  should  be  noted  that,  being  statistical  in  nature,  the  LPD  method  enables  the  refinement 
of  the  histogram  and  improvement  of  the  accuracy  by  subsequent  accumulation  of  bubble 
statistics.  Since  the  statistical  error  is  proportional  to  1  In'12 ,  where  n  is  the  number  of  bubbles  in 
a  histogram  box,  it  would  require  four  times  as  many  bubbles  for  a  two-dimensional  histogram  to 
double  the  resolution  along  each  axis,  or  to  reduce  the  statistical  error  by  half  with  the  same 
resolution.  The  choice  of  histogram  size  is  a  trade-off  between  the  spatial  resolution  and 
statistical  error. 

It  should  be  noted  that  the  statistical  nature  of  the  Lagrangian  approach  can  be  of  an 
especial  benefit  when  using  distributed  memory  computer  platforms  to  compute  multiphase 
flows.  By  simply  running  the  same  simulation  on  several  processing  nodes  with  a  different  initial 
conditions  one  can  increase  the  sample  size  and  improve  the  statistics.  A  particularly  attractive 
feature  of  this  approach  is  that  no  domain  decomposition  will  be  necessary  for  these 
computations  and  therefore  a  strictly  linear  scalability  can  be  achieved.  This  is  the  case  of  what 
is  sometimes  referred  to  as  an  embarrassingly  parallel  computation. 

5.4  Ship-wake  Simulations 
5.4.1  LES  Results 

Ship-wake  simulations  were  setup  to  correspond  to  the  un-propelled  Navy  5415  model 
[Smirnov  et  al.,  2001],  The  Reynolds  number  based  on  the  ship  width  is  on  the  order  of  2  •  1 05 . 
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Computational  grid  with  dimensions  164x98x92  was  used  to  represent  the  near- wake  region 
(Fig.  5.7)  of  one  ship-length,  L ,  in  the  axial  direction,  0.3  L  in  depth  and  0.61  in  the  span-wise 
direction.  The  simulations  were  performed  on  a  1GB  1GHz  DEC  Alpha  computer,  with  an 
average  computing  time  of  one  month  for  a  single  flow-through  time. 

Only  the  wake  flow  was  simulated,  with  the  inlet  plane  conditions  provided  by  prior 
RANS  results  carried  out  at  the  university  of  Iowa  [Ramaprian,  1981].  A  distribution  of  turbulent 
kinetic  energy  in  the  inlet  plane  is  shown  in  Fig.  5.8.  The  turbulent  kinetic  energy  data  were  used 
to  set  the  amplitudes  of  unsteady  subgrid-velocity  field  by  means  of  RFG  procedure  [Smirnov  et 
al.,  2001].  At  this  stage  RFG!  was  not  used  for  stochastic  tracking  of  the  bubbles  inside  the 
domain  as  was  done  in  the  mixing  layer  simulations.  The  free  surface  was  assumed  to  be  flat, 
invoking  the  low  Froude  number  approximation. 

Bubble  motion  was  governed  by  the  model  of  Sridhar  and  Katz  (1).  Bubbles  of  100//  m 
in  size  were  randomly  injected  at  the  inlet  plane  with  the  probability  distribution  proportional  to 
the  turbulent  kinetic  energy  level  (Fig.  5.8).  At  the  free  surface  the  life-time  of  the  bubbles  was 
set  to  zero  for  the  absence  of  more  exact  empirical  data.  The  integration  of  (Eqn.  5.1),  via  a 
second-order  Runge-Kutta  scheme  provides  the  new  velocity,  v  (?) ,  in  the  x,  direction  for  each 
particle  as  a  function  of  time. 

It  should  be  noted  that  the  range  of  bubble  sizes  used  to  calibrate  Eqn.  5.1, 
500  -  700//  m ,  is  larger  than  the  bubble  size  used  in  simulations.  Although  our  validation  on  the 
case  of  mixing  layer  may  indicate  that  approximation  (Eqn.  5.1)  may  still  be  valid  for  smaller 
bubble  sizes,  a  more  through  investigation  of  the  influence  of  bubble  size  may  be  appropriate. 

Figure  5.9  shows  an  instantaneous  picture  of  all  the  bubbles  in  the  wake  as  observed  from 
the  rear-bottom  comer  of  the  computational  domain  (ship  stem  is  not  shown).  It  can  be  seen  that 
the  bubbles  tend  to  cluster  in  preferential  concentration  regions  corresponding  to  high  vorticity. 
These  regions  are  dictated  by  the  dynamics  of  the  flow  in  the  near  wake.  Figure  5.10  provides  a 
cross-sectional  view  of  an  instantaneous  velocity  and  bubble  distributions  in  the  wake,  which  is 
rather  typical.  Two  symmetrical  vortex  structures  in  the  velocity  distribution  coincide  with  the 
bubble  agglomeration  regions.  With  the  higher  vortex  intensities,  which  will  take  place  at  a 
higher  Reynolds  number,  bubble  distribution  in  the  wake  will  be  increasingly  affected  by  these 
vortices.  This  is  especially  true  for  smaller  bubbles. 
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Figure  5.11  shows  typical  cross-sectional  distributions  of  bubbles  in  the  wake.  These 
distributions  are  not  instantaneous,  but  represent  total  counts  of  bubbles  accumulated  over  the 
simulation  run. 

Figure  5.12  shows  the  contours  of  probability  density  functions  of  bubble  occurrence  in 
the  wake.  These  probabilities  were  computed  from  bubble  distribution  histograms  obtained  in  the 
similar  way  as  for  the  mixing  layer  validation  case  described  above.  The  profiles  shown  in  Fig. 
12  indicate  a  rapid  bubble  population  decay  and  gradual  spreading  of  the  bubble  cloud  at  the 
distance  of  one  ship-length.  Although  the  classified  nature  of  bubble  measurement  data  in  the 
wakes  of  Navy  ships  prevents  us  from  making  a  direct  comparison,  we  found  that  the  predicted 
bubble  distributions  are  similar  to  those  observed  in  typical  ship  wakes  [Hyman,  1998;  Hyman, 
2000]. 

Table  5.1  shows  the  total  bubble  counts  in  different  planes  Ci=u  and  the  corresponding 
normalized  ratios:  R =  C, / C, .  These  distributions  indicate  a  strong  depletion  of  bubbles  over  the 

half  ship-length  distance,  shown  in  Fig.  5.13.  The  depletion  of  bubbles  is  due  to  the  buoyancy 
driven  migration  to  the  surface  and  the  dissolution  effects.  Buoyancy  forces  affect  mainly  large 
bubbles,  whereas  the  dissolution  affects  primarily  small  ones.  Therefore,  even  the  small  bubbles 
entrained  by  the  vortexes  will  eventually  disappear  from  the  domain. 


PLANE 

COUNT 

RATIO 

1 

46998 

0.01 

2 

9803 

0.21 

3 

3326 

0.07 

4 

1119 

0.02 

Table  5.1  Cumulative  bubble  distributions  in  different  planes 

5.4.2  Bubble  Clustering  Effects 

Bubbles  motion  in  a  turbulent  vortex  may  lead  to  clustering  in  the  vortex  center.  This,  in 
turn,  may  cause  bubble  coalescence  and  influence  the  flow  itself.  In  order  to  analyze  the 
conditions  under  which  such  clustering  takes  place  we  considered  the  case  of  a  realistic  vortex 
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flow  obtained  from  experimental  data  [Larsson  et  al.,  1998].  An  equation  of  bubble  motion, 
including  all  the  important  terms  [Sridhar  and  Katz,  1995]  was  used  to  predict  bubble  behavior. 
The  results  were  analyzed  to  determine  the  limiting  values  of  non-dimensional  parameters  under 
which  bubble  clustering  and  coalescence  effects  may  become  important. 

The  objective  of  the  computations  was  to  capture  the  regime  where  the  bubble  would 
converge  to  the  center  of  the  vortex  rapidly  enough  to  make  the  coalescence  events  more  likely. 
Then,  using  two  non-dimensional  parameters,  a ,  representing  the  vortex  strength  (Stokes 
number)  and  p ,  representing  the  ratio  of  bubble-size  to  bubble-separation,  general  conditions  for 
bubble  clustering  can  be  obtained. 

Given  two  bubbles  separated  by  a  distance  d  in  a  vortex  centered  at  the  middle  point 
between  the  bubbles,  the  necessary  condition  for  bubble  coalescence  is 

CR>P~'  (5.6) 

where  CR  is  a  convergence  rate,  defined  as 

c„=—  (5.7) 

ri 

with  r0  is  the  initial  distance  of  the  bubble  from  the  axis  and  rx  its  distance  from  the  axis 
after  one  revolution  in  a  vortex.  Condition  (Eqn.  5.6)  means  that  the  bubbles  should  come  into 

n 

contact  with  each  other  after  one  revolution  in  the  vortex  . 

Coalescence  criterion  (Eqn.  5.6)  may  be  rather  strong.  In  reality  bubble-bubble 
interaction  effects  may  become  important  long  before  the  criterion  (Eqn.  5.6)  is  satisfied. 
Consequently,  we  introduce  the  criterion  of  bubble  clustering 

CR>0.5p-'  (5.8) 

which  identifies  the  bubbles  that  approach  each-other  to  the  distance  of  one  bubble- 
diameter  or  less. 

To  compute  bubble  motion  in  the  turbulent  vortex  equation  (5.1)  was  used  to  describe  the 
bubble  motion  and  the  flow-field  was  constructed  using  a  curve-fit  to  the  experimental  data 
[Larsson  et  al.,  1998]  (Fig.  5.14).  The  following  function  was  used  to  approximate  the  data 
Uz=a(\- Bz  exp(-(r /RJ2))  (5.9) 


7 


Usually  turbulent  eddies  do  not  last  for  more  than  one  revolution 
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Ue  =arexp(-(r/Re)r) 

where 
a= 1.0 


(5.10) 


Bz=  0.08 
Rz=  0.7 
R0=O.22 

r  =  0.58  (5.11) 

The  advantage  of  using  the  analytical  formula  is  that  all  the  velocity  derivatives  required 
in  the  bubble-equation  of  motion  can  be  computed  analytically,  thereby  minimizing  the 
discretization  errors.  So,  the  spatial  derivatives  needed  for  calculating  the  vorticity,  a>  in  (5.1), 
can  be  easily  computed  from  the  expressions  (5.9),  (5.10),  which  gives  only  two  non-zero 
derivatives 

8UZ  _2BZ  ar  c-(rlR,)2 
dr  R] 

=a  (1  -y  {R^-r  rr)e-(r/R<>Y 
dr 

An  important  measure  of  the  strength  of  the  vortex  is  it’s  turnover  time.  For  the  whole 
core  region  this  time  is  almost  the  same,  since  the  core  region  can  be  approximated  by  the  solid 
body  rotation 

Ue=a-r  (5.12) 

where  a  is  the  strength  of  the  vortex,  which  is  related  to  the  turnover  time  of  the  core  as 
Te -2 na'\  Since  data  in  Fig.  14  are  non-dimensional,  we  used  a  as  a  parameter  to  define  the 
vortex  strength.  The  vortex  size  parameters  Rg,Re  were  kept  constant,  which  corresponds  to  the 
vortex  core  size  of  de  ^0.5 m .  We  did  not  scale  down  the  vortex  size  since  the  results  of  Fig.  5.14 
were  obtained  for  a  large  vortex  produced  in  the  stem  of  the  ship  and  it  is  not  clear  if  they  can  be 
scaled-down  to  small  eddy  sizes.  All  computations  were  applied  to  the  case  of  d<de  and 

d,<zd,  i.e.  the  bubbles  were  considered  to  be  inside  of  the  vortex  core  and  bubble  size  much 

b  e  3 

smaller  than  the  typical  eddy  size.  For  a  given  bubble  separation  d  the  results  would  apply  to  all 
the  eddy  sizes  greater  than  d ,  and  having  the  same  vortex  strength  a .  The  neglect  of  eddy  sizes 
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smaller  than  the  separation  length,  i.e.  de<d ,  may  lead  to  some  overestimate  of  clustering 
effects  for  highly  turbulent  flows. 

Figure  5.15  shows  a  typical  convergent  trajectory  of  a  bubble.  In  this  case  a  bubble  of 
rb  =  700nm  and  rb  =  1 .03  - 1 0“4  v  was  injected  into  a  core  of  a  vortex  with  a  =  1  Os'1  at  re  =1  cm 
from  the  axis.  The  convergence  factor  was  close  to  Cr=2.  The  corresponding  non-dimensional 
parameters  for  this  case  are:  «  =  1.64-10~4  and  /?  =  0.07,  which  lies  outside  of  coalescence 
limits  defined  by  (Eqns.  5.6  and  5.8). 

To  visualize  the  results  the  computed  particle  trajectories  were  plotted  in  a  coordinate 
system  moving  along  the  vortex  axis  with  the  velocity  equal  to  the  center-line  axial  velocity  of 
the  flow:  Uz(r  =  0) .  In  Fig.  5.15  gravity  acceleration  vector  g  is  directed  along  the  axis.  The 
main  contribution  to  the  bubble  dynamics  in  this  flow  came  from  the  added  mass  force  F0 , 
responsible  for  bubbles  attraction  toward  the  the  vortex  axis,  and  the  gravity  force,  FA ,  which 
caused  them  to  accelerate  along  the  axis.  As  can  be  seen  in  the  figure,  the  displacement  along  the 
axial  direction  is  smaller  than  the  drift  toward  the  center,  which  approximately  relates  to  the  ratio 
of  buoyancy  to  added-mass  forces  characteristic  for  this  particular  bubble  size  and  vOrtex 
strength.  In  only  about  1/3  of  all  cases  the  vortex  axis  can  be  considered  as  approximately 
aligned  with  the  gravity  force,  whereas  in  the  other  2/3  of  the  cases  the  corresponding  axes  are 
misaligned,  and  bubbles  convergence  toward  the  center  is  slower.  So,  for  the  bubble  with  CR  =  2 
and  axis  aligned  with  the  gravity  the  convergence  rate  went  down  to  approximately  1.5  when  the 
axis  was  turned  at  the  right  angle  to  g .  Considering  this,  a  more  accurate  convergence  criterion 


can  be  formulated  as 

(5.13) 

where 

CR  =  1/3  CR  ||  +  2/3  CR1 

(5.14) 

with  Cm  being  the  convergence  rate  for  the  parallel  and  CR1  for  the  normal  alignment  of 

the  gravity  and  vortex  axis.  As  can  be  seen,  large  enough  bubbles,  experiencing  strong  gravity, 
can  be  completely  drifted  away  from  the  vortex  in  the  normal  alignment.  It  should  be  noted  that 
the  coalescence  criteria  (Eqn.  5.6),  (Eqn.  5.8),  (Eqn.  5.13)  are  rather  approximate,  and  should 
only  be  used  for  order  of  magnitude  estimates. 
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Our  calculations  of  bubbles  with  db  =\mm  (rb  =  5.26-10~5)  show  (Fig.  5.16)  that  for 
P  -  0.25  the  coalescence  criterion  (Eqn  5.6)  is  approximately  satisfied  for  the  vortex  strength  of 
a  =  40C1  and  clustering  criterion  (Eqn  5.8)  -  for  a  =  20s~x .  This  corresponds  to  the  Stokes 
numbers  of  a> 3.34 •  1 0-4  and  a>1.67-10~4  respectively. 

n 

In  our  LES  computations  of  a  wake  behind  a  Navy  ship  model  the  vorticity  varies 
between  0  and  100.  Using  the  relation  \co\~a ,  valid  in  the  vortex  core  region,  we  can  get  an 
estimate  of  the  range  of  Stokes  numbers  that  may  be  encountered  in  simulations.  For  a  bubble  of 
db  =1  mm  it  is  0<a<8.33T0-4 .  For  bubbles  of  smaller  sizes  the  lower  bound  on  a  will  be 

lower.  Hence,  if  we  restrict  our  analysis  to  the  bubbles  of  sizes  less  than  1mm  in  size,  we  may 
infer  that  the  maximum  Stokes  number  that  we  can  currently  reproduce  in  computations  is 
a « 8.33  -10'4.  Considering  our  clustering  bound  for  bubble  coalescence  of  o:>1.67-10-4,  we 
may  conclude  that  in  our  current  LES  computations  we  may  consider  bubble  clustering  and 
coalescence  effects  in  the  range  of  (1.67<a<8.33)-10'4 . 

The  important  spatial  measure  is  the  distance  of  the  bubble  from  the  vortex  axis  rather 
than  the  size  of  the  eddy,  which  can  be  much  bigger.  Considering  this,  the  bubble/eddy  size  ratio, 
/?,  would  greatly  depend  on  bubble  size  distribution  and  concentration.  The  value  P-0.25 , 
which  was  used  in  the  previous  calculations,  requires  bubble  separation  from  the  axis  to  be  of 
4mm  or  less  for  all  bubble  sizes  under  1mm.  Consequently,  for  the  clustering  of  the  bubbles  to 
occur  their  concentration  should  be  high  enough  to  allow  at  least  two  bubbles  at  distances  shorter 
than  4mm  from  each  other,  which  gives  the  lower  bound  on  bubble  number  density  as 
Nb  >2/(4-l(T3)3  =  3.13 -107  m ~3 .  Combining  with  our  previous  estimate  of  a  ,  we  conclude  that 

for  bubble  number  densities  above  3.13 TO7  m~3  bubble  coalescence  effects  can  be  present  for 
a  >1.67 -10“4. 

On  the  other  hand,  considering  the  maximum  vorticity  in  the  wake  amax  =1005"' ,  we  may 

give  an  estimate  of  lower  bubble  concentration  limit  when  the  coalescence  effects  become 
important.  The  computations  showed  that  the  coalescence  criterion  (Eqn  5.6)  for  a  bubble  of 


8 
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db  =  1  mm  ( rb  =5.26-10  5)  in  a  vortex  core  of  a  =  1005  1  gives  >5  =  8.33  10  2 ,  which  provides 
the  following  coalescence  conditions  on  bubble  number  density,  Nb ,  and  void  fraction,  Fb : 

Nb°  > — — t=1.15-106  m"3 

WP) 

A  77" 

Fb°  > Nb°  —  r63  =  6.02  •  1 0“4  (5.15) 

Bubble  clustering  criterion  (Eqn  5.8)  results  in  /?  =  4.17  -10~2 ,  and  the  corresponding 
relations  for  the  number  density  and  the  void  fraction  are 

Nb  =  1/2 Nc°  >  — - — - = 5.75  •  1 05  m"3 

wpf 

Fbl  =  1/2 Fb°  >  3 .0 1 T  0”4  (5.16) 

With  the  higher  resolution  of  LES  than  we  are  currently  using,  and  which  we  shall  be 
able  to  ‘attained  on  distributed  memory  platforms,  higher  levels  of  vorticity  can  be  reached.  This 
is  usually  a  result  of  smaller  eddies  generating  higher  vorticities  than  parent-eddies.  It  may  bring 
the  estimates  (Eqn.  5.15),  (Eqn.  5.16)  further  down  in  number  densities  and  void  fractions. 

The  cases  above  relate  to  the  situation  of  mono-dispersed  bubbles  (single-size).  In  the 
poly-dispersed  case  the  coalescence  criteria  can  be  formulated  by  redefining  (5  as 

P=(db}+db2)/(2de)  (5.17) 

where  dbvdb2  are  diameters  of  bubbles  of  two  different  sizes  in  the  vortex,  separated  by 
de .  This  expression  for  ft  should  now  be  used  in  (Eqn  5.6)  and  (Eqn.  5.8).  To  obtain  an  estimate 
of  bubble  collision/coalescence  frequency  a  double  integration  over  all  the  bubble  sizes  dbvdb2 

weighted  by  the  bubble-size  distribution  would  be  required.  A  rough  estimate  can  still  be  made 
by  single  integration  of  the  results  for  the  mono-dispersed  case,  considering  that  the  obtained 
void  fraction  limits  will  represent  an  over-estimation  of  the  real  case. 

Since  estimates  (Eqns.  5.15  &  5.16)  were  based  on  LES  computations  of  a  ship  model 
with  a  relatively  small  length  of  5.72 m  we  can  expect  the  upper  bounds  on  vorticity  to  further 
increase  in  the  case  of  a  real  ship  where  the  Reynolds  number  and  vorticity  are  much  higher. 
Consequently,  the  estimates  of  Nb  and  Fb  in  (Eqns.  5.15  &  5.16)  can  be  lower. 
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We  should  also  note  that  clustering/coalescence  effects  considered  here  were  attributed  to 
a  single  revolution  of  the  vortex.  In  reality  an  accumulated  effect  of  many  revolutions  in  one  or 
many  subsequent  eddies  may  significantly  amplify  the  effect  of  a  single  revolution. 
Consequently  it  will  lead  to  lower  values  of  Nb,Fb  than  estimated  in  (Eqns.  5.15  &  5.16).  On  the 

other  hand  the  estimates  were  made  for  the  case  of  two  bubbles  aligned  on  the  opposite  side  of 
the  vortex  axis,  which  will  be  true  for  only  a  fraction  of  all  possible  bubble  alignments.  The 
account  of  this  will  increase  the  values  of  Nb,Fb. 

Considering  all  these  factors,  the  actual  values  of  void-fraction  limits  for 
clustering/coalescence  can  be  close  to  the  maximum  values  observed  in  real-size  experiments, 
which  is  about  10  4 ,  and  lie  within  the  range  of  values  that  can  be  expected  in  the  stem  region 
[Hyman,  2000].  The  numerical  studies  of  Carrica  et  al.  (1998)  provide  the  estimates  of  maximum 
void  fractions  of  1 0  3 ,  which  occur  mainly  near  the  free  surface  in  the  stern  region. 

Therefore,  we  can  conclude  that  in  realistic  ship-wake  situations  bubble  clustering  and 
coalescence  effects  may  be  important,  especially  in  the  near- wake  region  of  1-3  ship  lengths. 

5.5  Discussions 

The  results  of  this  study  show  the  viability  of  joint  LES/LPD/RFG  method  for  computing 
turbulent  bubbly  wakes,  which  can  be  applied  to  high  Reynolds-number  ship-wake  flows. 
Bubble  distributions  for  mixing  layers  and  wakes  were  obtained. 

The  effect  of  bubble  influence  on  the  flow  field  was  not  considered  in  this  study. 
Nevertheless,  the  work  by  others  [Elghobashi  and  Truesdell,  1993;  Truesdell,  and  Elghobashi, 
1994]  indicates  that  when  there  is  a  large  density  ratio  between  the  phases,  as  considered  in  this 
study,  the  influence  of  bubbles  on  the  carrier  phase  may  become  important  even  though  the  void 
fraction  remains  small.  Works  on  this  topic  are  continuing  under  the  current  project. 

To  further  improve  the  accuracy  of  predictions  for  the  flat-plate  wake  a  more  realistic  no¬ 
slip  wall  boundary  condition  may  be  used  together  with  the  grid  refinement  at  the  wall  so  as  to 
resolve  the  turbulent  boundary  layer. 

Stochastic  bubble  tracking  with  RFG  for  the  ship-wake  case  will  have  to  be  introduced 
and  compared  to  the  current  approach  and  experimental  data  to  build  a  more  thorough  ship-wake 
validation  case. 


106 


Performed  estimates  of  bubble  clustering  effects  show  that  they  can  be  important  in 
certain  regions  of  the  near  wake.  In  the  first  approximation  the  account  of  these  effects  can  be 
done  by  introducing  a  two  way  coupling  (bubble-flow  influence)  scheme.  In  the  second  level  of 
approximation  it  can  be  done  by  including  bubble  coalescence  effects,  which  can  be 
accomplished  within  an  implicit  probabilistic  interaction  scheme  proposed  earlier  by  the  authors 
[Smirnov  and  Celik,  2000], 

Given  long  execution  times  on  a  single  work  station,  it  would  be  appropriate  to  consider 
an  "embarrassingly  parallel"  run  on  a  a  computer  cluster  (Beowulf),  where  each  work  station 
would  run  essentially  the  same  simulation  with  the  initial  conditions  taken  from  independent 
statistical  samples.  This  way  a  linear  speedup  can  be  achieved  in  collecting  the  statistics  and 
improving  the  accuracy  of  the  computations.  This  assertion  has  been  investigated  and 
encouraging  results  have  been  obtained  for  shipwake  simulations  (see  Appendix  A). 
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Figure  5.3  Instantaneous  bubble  distribution 
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Figure  5.6  Bubble  concentrations  (^=0.079365  m). 
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THE  SKETCH  OF  FLOW  PAST  SHIP  HULL 
(L  is  the  ship  length) 

Figure  5.7  Near  and  far  ship  Wake  regions. 


Figure  5.8  Turbulent  kinetic  energy  in  the  coriiputational  inlet  plane  (m2/s2). 


Ill 


Figure  5.9  Bubbles  in  a  ship-wake.  Rear  bottom  view.  The  upper  surface  of  the  box  represents  the 
free  surface.  The  middle  rectangle  represents  the  inlet  plane.  More  detailed  cross-sectional  distributions 

are  shown  in  Figure  1 1 . 
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(b)  Bubble  distribution 


Figure  5.10  Velocity  vectors  and  bubble  distribution  at  X/L  =  0.2. 


Figure  5.15  Trajectory  of  a  bubble  in  a  turbulent  vortex  (R=700p.m). 
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Figure  5.16  Convergent  bubbles  (Bubble  motion  shown  is  relative  to  an  observer  moving  with  the 

axial  velocity  of  the  vortex.) 
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6  PARALLELIZATION  OF  THE  LES  CODE  FOR  SHIP  WAKES 


Fast  computers  have  stimulated  the  rapid  growth  of  a  new  way  of  doing  science.  The  two 
broad  branches  of  theoretical  science  and  experimental  science  have  been  joined  by 
computational  science.  Computational  scientists  simulate  on  supercomputers  phenomena  too 
complex  to  be  reliably  predicted  by  theory  and  too  dangerous  or  expensive  to  be  reproduced  in 
the  laboratory.  There  is  always  a  demand  for  greater  computational  speed.  Areas  requiring  great 
computational  speed  include  numerical  modeling  and  simulation  of  scientific  and  engineering 
problems.  The  cost  of  an  advanced  single  processor  computers  increases  more  rapidly  than  their 
power,  which  made  parallel  processing  evolve. 

In  this  part  of  the  study  we  present  the  parallel  implementation  of  large  eddy  simulations 
(LES)  of  a  ship  wake  using  the  domain  decomposition  technique.  We  present  the  results  of  the 
implementation  executed  on  a  cluster  of  workstations  at  West  Virginia  University  and  Pittsburgh 
Supercomputing  Center.  Also,  we  will  show  how  the  implementation  scales  up  with  the  number 
of  workstations  and  that  it  is  possible  to  obtain  better  accuracy  by  increasing  the  number  of 
workstations  in  the  cluster  system. 

Algorithmic  improvements  and  faster  machines,  particularly  parallel  machines,  provide 
the  opportunity  for  effectively  using  large-eddy  simulations  for  the  problems  of  practical 
importance.  The  main  objective  in  this  study  is  to  predict  the  development  of  turbulence  in  a  ship 
wake  flow  using  cluster  computing. 

This  requires  analysis  of  the  ship  wake  flow,  development  of  an  efficient  and  accurate 
simulation  of  turbulence  in  the  ship  wakes  by  refining  the  large  eddy  simulation  methodology, 
setting  up  the  computational  domain,  analysis  of  the  computer  resources,  and  analysis  of  the 
results. 

6.1  Methodology 

The  conventional  domain  decomposition  technique  for  elliptic  problems  is  realized 
through  a  two-way  exchange  of  data  at  the  boundaries  of  the  domains  (Simon,  1992;  Dihn  et  al., 
1984)  as  illustrated  in  Fig.  6.1(a).  This  guarantees  the  convergence  to  the  corresponding  single 
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domain  case.  However,  this  strategy  may  carry  an  excessive  communication  overhead  for  three 
dimensional  CFD  simulations.  If  the  problem  is  parabolic  ill  one  of  the  spatial  directions  one  can 
employ  a  parabolic  communication  approach  (Fig.  6.1(b));  This  may  reduce  communication 
overhead  by  half.  To  test  the  validity  of  this  approach  a  parallel  version  of  the  LES  code  has 
been  implemented  using  a  one-way  data  exchange  (Osman  dt  al.,  2000). 


(a)  Elliptic  decomposition 


(b)  Parabolic  decomposition 


(c)  Boundary  conditions 


Figure  6.1  Domain  decomposition  Strategy. 


The  drawback  of  the  parabolic  exchange  scheme  is  the  necessity  to  provide  additional 
outlet  boundary  conditions  for  each  domain,  which  can  altef  the  character  of  the  flow  close  to  the 
domain  outlet.  To  avoid  the  influence  of  this  distortion  oh  the  flow-  field  the  communication 
plane  should  be  set  at  some  distance  from  the  outlet  plane.  Thus  some  of  the  memory  space  and 
the  processing  time  is  inevitably  lost.  Moreover,  only  by  using  elliptic  message  transfer  can  one 
apply  domain  decomposition  technique  in  non-paraboliC  directions  which  is  necessary  for 
geometrically  complex  flows. 


120 


For  a  linear  solver  a  two-node  overlap  communication  strategy  (Fig.  6.1a)  would  be 
enough.  For  a  high-order  accurate  solver,  however,  a  wider  overlap  in  communication  planes 
will  be  appropriate.  Since  the  LES  solver  is  usually  at  least  second  order  accurate,  we  employed 
four  and  six-node  overlaps  in  our  communication  schemes.  For  second  order  accurate  diffusion 
terms  in  momentum  equations  a  four-point  overlap  is  enough,  whereas  for  pressure  solver 
decomposition  a  six  point  overlap  would  be  more  appropriate. 

6.2  Applications 

The  parallelizing  methodology  described  above  was  implemented  on  a  Beowulf  cluster  at 
WVU  and  at  Pittsburgh  supercomputer  center  (www.psc.edu).  Both  clusters  run  Linux  operating 
system  with  the  interprocessor  communications  based  on  MPI. 

6.2.1  One  Way  Communication  Simulations 

In  order  to  test  the  viability  of  the  implemented  decomposition  scheme  several  test 
simulations  were  performed.  The  communicated  data  were  velocities  and  contravariant  velocities 
and  pressure.  It  should  be  noted  that  since  the  pressure  solver  is  usually  sub-cycled  to  the  flow 
solver  sending  the  pressure  information  carries  most  of  the  communication  overhead.  In  practice, 
it  appears  that  excluding  pressure  communication  does  not  lead  to  a  significant  degradation  of 
the  solution  for  non-pressure  driven  flows,  like  wake  flows.  This  conclusion  was  also  confirmed 
in  our  simulations.  A  flat  plate  wake  flow  was  used  as  the  first  case.  The  geometry  and  numerical 
scheme  can  be  found  in  (Smirnov  et  al.  2001).  For  the  sake  of  a  higher  accuracy  the  simulations 
were  continued  with  two-way  communications  as  required  in  elliptic  solutions. 

6.2.2  Two  Way  Communication  Simulations 

Several  medium  scale  simulations  were  performed.  First  a  laminar  channel  flow  case  was 
computed  on  2,  4  and  8  processors.  Figure  6.2  shows  the  axial  velocity  contours  obtained  from 
the  8-processor  run.  Although  pressure  communication  was  blocked  in  this  simulation  the 
contour  lines  show  no  irregularities  at  the  inter-processor  boundaries.  These  results  confirm  the 
correctness  of  the  decomposition  scheme,  and  support  our  assumption  that  pressure  coupling  is 
rather  week  between  the  processors.  For  the  case  of  a  wake  flow  the  effect  of  pressure  will  be 
even  smaller,  thus  justifying  the  velocity-only  decomposition  strategy.  Despite  this  finding,  for 
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Figure  6.3  Comparison  of  numerical  meshes  used  in  single  arid  3-processor  simulation  (black  lines 

indicate  inter-processor  boundaries.) 

The  flat-plate  shear  layer  was  investigated  for  verification.  The  mesh  created  for  this  task 
is  shown  in  Figure  6.3.  Keeping  the  domain  size  and  effective  number  of  grid  points  (570K 
nodes)  constant,  the  same  simulation  was  performed  on  a  sihgle  processor  and  on  3  processors. 
The  comparison  of  the  predictions  as  shown  in  the  streamwise  contours  in  Figure  6.4  and  profile 
along  the  x-axis  in  Figure  6.5  has  shown  that  both  simulations  are  in  perfect  agreement.  This  has 
ensured  us  with  confidence  in  the  predictions  of  parallel  sithulations  of  this  type  of  flows.  The 
processor  time  was  also  recorded  for  this  case  on  the  local  DecAlpha  Cluster  (570K  grid  nodes). 
Here,  the  single  processor  execution  took  7.15  seconds  per  pits  (pits  =  Pressure  Iteration  *  Time 


Step)  and  the  3  processor  execution  took  2.78  sec/pits  per  CPU,  which  indicated  a  scale  up  of 
2.57  or  85%  efficiency. 


Figure  6.4  Comparison  of  predicted  streamwise  velocity  contour  in  single  and  3-processor  simulations 

(black  lines  indicate  inter-processor  boundaries.) 


Figure  6.5  Comparison  of  predicted  streamwise  velocity  at  one  point  in  single  and  3-processor 

simulations 

Next,  two  flat-plate  wake  simulations  on  4  and  8  processors  were  done  for  the  wake  flow 
of  Reynolds  number  1.2-106.  In  both  simulations  the  total  number  of  grid  nodes  was  equal  to 
224x18x10,  with  28x18x10  nodes  per  processor  in  8-processor  run  and  56x18x10  nodes  per 
processor  in  a  4-processor  run.  The  results  were  compared  to  look  for  any  possible  discrepancy 
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introduced  by  the  inter-processor  communications.  Figure  6.6  shows  a  perfect  agreement 
between  these  two  cases  on  the  computed  flow-field.  The  virtual  absence  of  turbulent 
unsteadiness  in  the  figures  is  the  result  of  high  numerical  diffusion  related  to  relative  coarseness 
of  the  grid.  Another  simulation  of  a  shear  layer  flow  was  performed  on  8  processors  with  the  grid 
size  of  250K  nodes  on  each  processor.  The  maximum  Reynolds  number,  based  on  shear  layer 
thickness  was  375.  As  can  be  seen  in  Fig.  6.7  the  development  of  shear-layer  was  not  affected  by 
inter-processor  communications.  The  speed  of  execution  of  the  large-scale  run  was  considerably 
slower  with  one  iteration  computed  in  approximately  1 .5  sec  on  a  DEC-Alpha  cluster  consisting 
of  533MHz,  512MB  nodes. 

The  scalability  analysis  performed  for  different  domain  decompositions  (Osman  et  al., 
2000)  indicted  that  the  speedup  is  almost  linearly  proportional  to  the  number  of  processors 
(domains)  being  used. 


Figure  6.6:  Comparison  of  wake  simulations  on  different  number  of  processors. 
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Figure  6.7:  Large-scale  8-processor  simulation  (25 OK  nodes  on  each  processor). 


6.2.3  Parallel  LES  Simulations  of  a  Ship  Wake  on  a  Straight  Track 

The  case  presented  in  section  4.3  is  considered  for  the  parallel  simulations  in  this  part  of 
the  study.  In  order  to  avoid  solving  for  the  flow  around  the  ship  model,  the  computational 
domain  starts  from  the  inflow  boundary  (or  initial  data  plane)  located  immediately  after  the  body 
in  the  wake,  where  there  is  no  flow-reversal.  The  pseudo-fandom  flow  field  generated  by  the 
RFG  method  is  added  to  the  mean  flow  of  RANS  simulatidhS  to  establish  the  unsteady  boundary 
conditions  at  the  inlet  plane.  The  whole  ship  flow  including  the  ship  wake  can  be  sketched  as  in 
Fig.  4.1(b).  The  computational  domain  starts  from  x/L  =  1:05  (where  x  starts  from  the  front  of 
the  ship  model).  At  this  plane  the  RFG  method  is  used  in  conjunction  with  the  RANS 
calculations  (Stem  and  Wilson,  2000). 

At  the  inflow  boundary  all  components  of  the  velocity  are  specified  as  a  function  of  time 
and  space.  At  the  outflow  boundary  Neumann  (free  gradiiflt)  boundary  conditions  are  applied. 
Symmetry  conditions  have  been  used  in  y  direction  and  periodic  boundary  conditions  have  been 
used  in  the  spanwise  (z)  direction.  At  the  free  surface  a  slip  condition  is  allowed  in  x  and  z 
directions,  but  the  velocity  component  normal  to  the  free  surface  is  set  to  zero.  As  such  the  free 
surface  is  approximated  as  a  moving  flat  plane,  i.e.  loW  Froude  number  approximation  is 
involved. 

The  computational  domain  size  is  3.0x0.3x0.5  (given  in  non-dimensional  units  in  ship 
length)  in  x,  y  and  z-directions  (axial,  vertical  and  transverse  directions),  respectively.  The  grid 
size  is  108x50x66  per  processor  multiplied  by  10  processors,  which  sums  up  to  ~3  million  grid 
nodes.  Non-uniform  grid  spacing,  stretching  smaller  tMti  1:03,  is  used  in  both  z-  and  y- 
directions.  The  length  scale  and  time  scale  used  in  RFG  wire  selected  as  constant  in  this  case. 
The  length  scale  was  0.02  of  the  ship  length,  and  the  time  siale  was  0.001,  non-dimensionalized 
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by  free  stream  velocity  and  ship  length.  Those  numbers  were  selected  because  the  turbulence 
length  scale  is  estimated  to  be  about  15%  of  the  ship  width.  The  Smagorinsky  SGS  model,  and 
second  order  central  differencing  scheme  were  used  unless  otherwise  stated. 

The  unsteady  streamwise  velocity  variation  is  presented  in  Figure  6.8,  where  x  is  the  flow 
trough  time.  This  point  was  selected  inside  the  wake,  as  seen  from  the  velocity  defect.  Figures 
6.9  and  6.10  show  the  vertical  variations  of  the  mean  velocities  and  the  root  mean  square 
velocities,  respectively,  where  the  time  averaging  was  made  over  3x.  Figure  6.9  shows  that  the 
spanwise  velocity  at  the  center  does  not  change  much  in  the  axial  direction,  but  same  velocity 
components  off  centerline  become  more  uniform  as  the  wake  develops.  The  vertical  velocity  at 
the  centerline  has  a  peak  below  the  free  surface,  which  gets  weaker  in  the  axial  direction. 
Significant  changes  are  observed  in  the  streamwise  velocity  profiles  in  both  axial  and  vertical 
direction.  The  location  of  peak  axial  velocity  moves  deeper  into  the  wake  indicating  a  plunging 
effect  as  the  wake  develops  in  the  streamwise  direction.  This  peak  is  stronger  off  the  centerline 
indicating  the  strength  of  one  of  the  bilge  vorticities  seen  in  Fig.  6.1 1.  The  intensity  of  the  rms 
velocity  fluctuations  show  very  little  decay  but  significant  redistribution  and  become  more 
isotropic  as  the  wake  develops.  Initially  the  turbulence  is  concentrated  more  near  the  free 
surface.  This  area  plunges  deeper  with  increasing  axial  distance. 

The  development  of  the  wake  in  the  axial  direction  is  presented  in  Figures  6.1 1  and  6.12. 
The  two  large  bilge  vortices  seem  to  first  move  away  from  each  other,  then  further  downstream 
they  gradually  merge  together  as  they  loose  strength.  However,  the  resolved  turbulent  kinetic 
energy  seems  to  increase  in  certain  regions  as  seen  in  Figures  6.1  lb  and  6.12b,  even  towards  the 
end  of  the  domain  (see  also  Fig.  6.15).  This  indicates  that  the  domain  should  be  longer  than  the 
selected  3  ship  lengths.  Also,  one  has  to  account  for  the  wake  spreading  and  a  deeper  and  wider 
cross-section  is  necessary  to  accommodate  the  whole  wake  further  downstream  than  3  Ship 
lengths. 

In  addition  to  these,  the  wake  spreading  or  wake  width,  which  is  obtained  from  the 
parallel  run  to  be  roughly  w~x1/4  is  also  consistent  with  Buller  and  Tunaley  (1989)’s 
measurements.  Milgram  et  al.  (1993)  and  Hoekstra&  Ligtelijn  (1991)  found  w~x1/5.  And  the 
spreading  rate  of  the  ship  wake  on  a  straight  track  is  observed  to  be  significantly  higher  than  that 
of  a  ship  wake  on  a  circular  track.  This  is  well  observed  from  the  parallel  simulation  predictions 
of  the  time  averaged  streamwise  velocity  magnitude  contours  in  Figure  6.11  for  a  straight  track 
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ship  wake  using  3  million  grid  nodes  on  10  processors  (Intel  Pentium4)  on  a  domain  of  3  ship 
lengths.  These  parallel  computations  were  also  expanded  to  a  6  million  nodes  simulation  on  6 
processors  (3  rd  generation  Dec  Alpha),  and  more  detailed  structures  were  captured.  However,  the 
turn  around  time  for  these  simulations  has  significantly  increased.  The  predicted  extra  detail  is 
especially  pronounced  in  the  vorticity  contours  presented  in  Figure  6.13.  The  6  million  nodes 
simulation  was  still  in  the  early  stages,  therefore  one  should  not  fall  into  the  wrong  conclusion, 
that  the  wake  spreading  is  much  more  pronounced.  In  fact,  it  will  tighten  as  the  simulations 
continue. 

The  resolved  turbulence  kinetic  energies  of  the  wake  simulations  (on  a  single  processor) 
of  a  ship  cruising  on  a  straight  and  circular  track  are  presented  in  Figure  6.14,  where  a  significant 
decay  of  kinetic  energy  in  the  streamwise  direction  can  be  observed.  The  locations  for  both 
studies  are  taken  to  be  at  the  highest  kinetic  energy  value  obtained  from  the  IDP  (along  the 
centerline  for  the  non  turning  ship  and  at  y=  -0.001  and  z=-3.18  for  the  turning  ship  wake 
simulation).  Comparing  to  the  predictions  of  a  non-turning  ship  wake,  the  turning  ship  case 
indicates  less  kinetic  energy  values.  However,  it  should  be  noted  that,  this  may  be  due  to  the 
coarser  grid  resolution  in  the  far  wake.  Overall  the  trends  are  similar,  but  there  is  a  sinusoidal- 
like  distribution  of  the  TKE  prediction  in  the  near  wake  for  the  straight  ship  case.  It  may  be 
because  of  the  existing  surface  wave  from  the  RANS  calculations  (Stem  and  Wilson,  2000).  This 
indicates  that  some  wave  information  may  be  present  implicitly  in  the  inflow  boundary. 
However,  this  sinusoidal-like  distribution  of  TKE  is  not  seen  in  the  simulations  of  a  ship  on  a 
circular  track,  which  may  imply  that  there  is  not  any  wave  information  present  in  the  RANS 
simulations  (Hyman,  2001). 


Figure  6.8.  Temporal  history  of  streamwise  velocity  components  at  x  =  0.6,  y  =  -0.012,  z  =  -0.06 
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Figure  6.9.  The  root  mean  square  velocity  profiles  at  3  spanwise  locations  at  a)  x’/L=0.09,  b)  x’/L=l  .60, 


c)  x’/L=2.80  (parallel) 


(a) 
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Figure  6.12  Contours  on  horizontal  y-planes  of  time  averaged  a)  streamwise  velocity  b)  turbulent  kinetic 
energy  (parallel  simulations,  108x50x66  nodes  per  processor) 
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Figure  6.13  a)  Vorticity  magnitude  contours  obtained  for  the  wake  behind  a  ship  on  a  straight  track  using 
parallel  computations  with  3  Million  grid  nodes,  b)  Preliminary  Vdrticity  magnitude  contours  obtained  for 
the  wake  behind  a  ship  on  a  straight  track  using  parallel  computations  with  6  Million  grid  nodes  (black 
lines  show  processor  boundaries,  overlap/communication  regions) 


Figure  6.14  Resolved  turbulence  kinetic  energy  of  the  wake  simulation  of  a  ship  cruising  on  a  straight  and 
circular  track  normalized  w.r.t.  its  inlet  value  (siftfle  node  computation) 


131 


Figure  6.15  Resolved  turbulence  kinetic  energy  of  the  wake  simulation  of  a  ship  cruising  on  a  straight 
track  normalized  w.r.t.  its  inlet  value  along  a  line  of  center  z  =  0.055,  y=-0.012  (Parallel  computations) 

The  same  dip  seen  in  Figure  6.14  is  also  present  in  the  3  million  nodes  simulation 
predictions.  However,  the  kinetic  energy  presented  in  Figure  6.15  also  indicates  that  the  kinetic 
energy  is  increasing  in  the  downstream  direction.  However,  the  increase  in  TKE  towards  the  end 
of  the  wake  is  an  art  effect  of  taking  a  line,  along  which  TKE  is  computed,  that  is  not  aligned 
with  the  main  streamwise  direction.  To  have  a  better  perspective,  Figures  6.1  lb  and  6.12b  should 
be  investigated.  Here,  the  bilge  vortices  separate  and  later  merge  closer  as  they  become  larger, 
which  is  probably  the  cause  for  the  turbulent  kinetic  energy  rise.  Moreover,  if  the  ship  turns  on  a 
circular  track,  the  wake  generated  by  the  ship  may  not  follow  a  circular  track.  Hence,  a  line  that 
is  not  aligned  with  the  center  of  a  bilge  vortex  would  give  a  false  idea  of  the  TKE  variation.  It 
should  be  mentioned  that  some  of  the  initial  decay  in  k  could  be  the  out  effect  of  the  application 
of  RFG  at  IDP.  This  issue  is  currently  under  investigation. 

6.3  Discussions 

Detailed  LES  calculations  have  been  performed  for  the  developing  wake  of  a  surface  ship 
cruising  on  a  straight  track  with  3  million  nodes  (108x50x66  nodes  per  domain,  Ax~0.46m; 
Ay~0.5m;  Az~1.2m  in  the  wake  region)  and  6  million  nodes  (330x50x66  nodes  per  domain, 
Ax~0.23m;  Ay~0.5m;  Az~1.2m  in  the  wake  region). 

The  results  are  analyzed  and  turbulence  statistics  have  been  presented  for  further 
comparison  with  RANS  and/or  experiments  as  they  become  available.  Qualitative  comparison 
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with  experimental  observations  indicate  that  the  LES  results  are  credible,  but  rigorous  validation 
is  necessary  to  make  definite  conclusions. 

Relatively  under-resolved  LES  of  the  wake  of  a  ship  turning  on  a  circular  track  has  also 
been  performed  and  the  results  are  compared  to  those  of  a  ship  cruising  on  straight  track. 

It  seems  that  the  wake  of  the  turning  ship  is  much  narrower  and  has  more  concentrated 
vorticity  due  to  the  merger  of  the  two  bilge  vortices.  There  is  little  similarity  between  these  two 
wakes.  The  wake  of  the  turning  ship  exhibits  more  dynamic  futures.  A  thorough  comparison  of 
the  two  cases  will  only  be  possible  when  the  grid  resolution  of  the  turning  ship  wake  is  increased 
by  employing  parallel  runs. 


133 


7  CONCLUSION 


Large  eddy  simulation  of  complex  turbulent  flows  and  bubble  dynamics  in  the  wake  of 
cruising  ships  has  been  studied  in  detail.  To  reduce  cost  and  simulation  time  the  calculations 
were  started  using  an  inflow  plane  in  the  wake  of  the  ship  itself.  This  requires  a  robust  approach 
for  providing  an  instantaneous  velocity  field  in  conjunction  with  the  prescribed  mean  flow  field. 

A  new  technique  was  developed  for  this  purpose,  namely  the  RFG  technique.  RANS  calculations 
which  were  including  the  ship  hull  were  used  to  provide  the  RFG  procedure  with  the  information 
needed  on  the  inflow  boundary.  The  further  development  of  the  wake  flow  was  calculated  via 
LES  using  the  domain  decomposition  technique  on  parallel  machines.  Furthermore,  the  flow  was 
injected  with  bubbles  to  compute  bubble  statistics 

Detailed  LES  calculations  have  been  performed  for  the  developing  wake  of  a  surface  ship 
cruising  on  a  straight  track  with  3  million  nodes  (108x50x66  nodes  per  domain,  Ax~0.46m; 
Ay~0.5m;  Az~1.2m  in  the  wake  region)  and  6  million  nodes  (330x50x66  nodes  per  domain, 
Ax~0.23m;  Ay~0.5m;  Az~1.2m  in  the  wake  region. 

The  results  are  analyzed  and  turbulence  statistics  have  been  presented  for  further 
comparison  with  RANS  and/or  experiments  as  they  become  available.  Qualitative  comparison 
with  experimental  observations  indicate  that  the  LES  results  are  credible,  but  rigorous  validation 
is  necessary  to  make  definite  conclusions. 

The  Random  Flow  Generation  (RFG)  technique,  that  has  been  developed  can  be  used  to 
initialize  the  turbulent  flow  field  and  also  to  prescribe  realistic  turbulent  inflow  boundary 
conditions  for  LES.  By  using  this  technique,  the  turbulent  flow  field  can  be  reproduced  with 
good  accuracy,  even  when  limited  information  is  available  on  the  mean  flow.  This  information 
can  include  the  turbulent  kinetic  energy,  its  dissipation  rate,  or  shear  stresses. 

The  wake  simulations  for  the  ship  model  DTMB  5512  moving  on  a  straight  track  has 
been  compared  with  the  macro  wake  measurements  by  Hoekstra  &  Ligtelijn  (1991).  The 
prominent  flow  structures  such  as  the  bilge  vortices  are  captured  by  LES.  In  measurements  and 
simulations,  the  minimum  axial  velocity  occurs  near  the  free  surface  of  the  center  of  the  wake. 
Moreover,  the  extent  of  axial  turbulence  intensities  is  in  good  agreement  with  measurements. 
Additional  smaller  side  vortex  pairs  are  observed  away  from  the  center  of  the  wake.  As  the  wake 
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widens,  the  strength  of  these  vortex  pairs  weaken.  The  vorticity  becomes  more  concentrated  near 
the  free  surface.  Moreover,  the  wake  spreading  or  wake  width,  which  is  obtained  to  be  roughly 
w~x1/4  is  consistent  with  Tunaley  &  Buller  (1989)’s  measurements. 

It  seems  that  the  wake  of  the  turning  ship  is  much  narrower  and  has  more  concentrated 
vorticity  due  to  the  merger  of  the  two  bilge  vortices.  There  is  little  similarity  between  these  two 
wakes.  The  wake  of  the  turning  ship  exhibits  more  dynamic  futures.  A  thorough  comparison  of 
the  two  cases  will  be  possible  when  the  grid  resolution  of  the  turning  ship  wake  is  increased  by 
employing  parallel  runs.  This  is  the  topic  of  an  upcoming  publication  which  will  be  an  addendum 
to  this  report. 

The  classical  SGS  model,  namely  the  standard  Smagorinsky  model,  is  not  totally  suitable 
for  complex  flows  as  it  uses  a  constant  eddy  viscosity  coefficient  for  the  entire  domain.  It  was 
observed  that  the  resolved  eddies  and  the  kinetic  energy  are  sensitive  to  the  eddy  viscosity 
coefficient.  To  remedy  this,  the  Smagorinsky  model  was  modified  to  account  directly  for  the 
effect  of  the  free  surface  on  turbulence  generation,  and  the  performance  a  nonlinear  one-equation 
SGS  model  was  assessed. 

It  has  been  clearly  demonstrated  that  if  there  are  wall  boundaries,  there  should  be  wall 
modifications.  The  results  indicated  that  very  reasonable  turbulence  statistics  were  predicted  with 
the  mentioned  modifications  and  it  could  be  used  as  a  SGS  model  in  the  simulation  of  ship 
wakes.  These  models  were  applied  to  the  wake  behind  a  circular  track  and  has  been  compared 
with  the  standard  Smagorinsky  model. 

The  results  of  this  study  show  the  viability  of  a  joint  parallel  LES/LPD/RFG  method  for 
computing  the  developing  turbulent  bubbly  wakes  generated  by  surface  ships,  at  high  Reynolds- 
numbers.  This  is  well  documented  in  our  publications  (see  References)  by  Celik,  Smirnov, 
Yavuz,  Shi,  Cehreli  and  Hu. 
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A  major  drawback  in  the  Lagrangian  particle  simulation  in  dispersed  two-phase  flows,  in 
terms  of  the  computational  cost  and  machine  capacity,  is  the  limitation  on  the  number  of 
particles,  or  particle  clouds  whose  trajectories  are  to  be  tracked  parallel  to  the  solution  of  the 
continuous  flow  field.  Insufficient  number  of  particles  being  tracked  commonly  leads  to 
inaccurate  statistics. 

Earlier  work  of  the  authors  on  parallelization  of  a  LES  solver  by  means  of  domain 
decomposition  [1,  2]  provided  the  possibility  to  simulate  large  scale  turbulent  structures  of 
typical  ship-wakes  on  computer  clusters.  The  next  logical  step  is  to  extend  the  pure  turbulence 
model  with  important  multi-phase  features  of  the  wake  such  as  bubble  dynamics.  The  algorithm 
for  particle  tracking  and  population  dynamics  developed  earlier  by  the  authors  demonstrated  the 
ability  to  efficiently  simulate  large  populations  of  particles  including  coalescence  effects  with 
even  modest  computer  resources  [3,  4,  5,  6].  However,  parallel  implementation  of  a  discrete 
particle  dynamics  algorithm  and  LES  flow  solver  by  means  of  domain  decomposition  technique 
commonly  leads  to  large  communication  overheads  and  load  balancing  problems.  In  this  study 
we  pursued  a  simple  embarrassingly  parallel  strategy,  which  enabled  us  to  avoid  these  two 
problems.  The  basic  idea  is  to  perform  the  simulation  of  statistically  independent  realizations  of 
the  flow-field  and  particle  ensembles,  with  each  realization  assigned  to  one  cluster  node.  This 
strategy  completely  excludes  any  communication  between  the  computing  nodes,  at  the  same  time 
achieving  the  perfect  load  balance.  However,  the  technique  calls  for  a  new  compromise  between 
the  required  accuracy  in  the  resolution  of  flow  features  and  the  desired  quality  of  particle 
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statistics.  With  the  whole  computational  domain  residing  on  one  node,  the  accuracy  of  flow 
solution  is  restricted  by  the  single  node  memory  size. 

Despite  its  seeming  simplicity  this  parallelization  strategy  should  take  care  of  ensuring 
statistical  independence  of  ensembles  generated  on  different  computing  nodes.  This  is  achieved 
by  imposing  independent  random  generation  of  three  important  flow  and  particle  conditions:  (1) 
initial  conditions  on  the  flow-field,  (2)  inlet  flow  conditions,  (3)  particle  injection  distribution. 
All  three  conditions  are  subjected  to  randomized  time-dependent  change,  which  nevertheless 
follows  a  predefined  statistical  distribution. 

It  should  be  noted,  that  while  generating  independent  ensemble  of  particle  injection 
coordinates/velocities  with  a  given  distribution  is  a  simple  matter  of  Monte-Carlo  sampling,  the 
generation  of  randomized  inflow  and  initial  conditions  for  the  flow  field  should  be  subjected  to 
certain  restrictions  imposed  by  continuum  dynamics  laws.  For  example,  the  continuity  relation 
will  generally  not  be  satisfied  for  any  random  distribution  of  flow  velocity  vectors,  even  if  this 
distribution  obeys  Gaussian  or  other  valid  statistics.  This  problem  of  adequately  representing 
time  varying  random  fluid  velocity  field  was  solved  earlier  by  the  authors  in  the  development  of 
RFG  method  [7],  and  was  applied  successfully  in  this  case  to  generate  statistically  independent 
and  divergencefree  flow-field  ensembles. 

Iterations  of  the  discrete  bubble  solver  were  sub-cycled  inside  flow  iterations  following 
Lagrangian  particle  dynamics  (LPD)  in  a  dilute  dispersed  two-phase  flow  with  one  way 
coupling,  i.e.  where  the  particles  motion  is  determined  by  the  continuous  phase  flow,  but  the 
continuous  flow  is  not  influenced  by  the  particles.  The  number  of  particles  tracked  in  each 
simulation  can  vary,  depending  on  the  machine  capacity  and  optimal  computational  expense.  The 
realizations  from  different  runs  on  different  nodes  were  collected  and  analyzed  to  produce 
histograms  of  bubble  distribution.  The  global  particle  statistics  was  obtained  by  averaging  over 
all  the  ensembles. 

The  combined  LES/LPD  solver  was  set  up  to  run  on  a  Beowulf  cluster,  with  1GHz  1GB 
computing  nodes.  A  set  of  simulations  of  a  turbulent  bubbly  ship  wake  flow  was  performed,  in 
which  a  total  of  about  254000  particles  were  tracked  on  nine  different  computing  nodes.  The 
combined  statistics  is  compared  to  the  statistics  from  a  single  run,  indicating  qualitatively 
equivalent,  but  much  better  results. 
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Figure  1  shows  a  snapshot  of  a  particle  ensemble  representing  bubbles  in  a  ship  wake. 
The  cumulative  cross-sectional  distribution  of  bubbles  in  the  inlet  plane  is  shown  in  Fig.2.  Only 
the  near-wake  regionwas  modeled  in  this  simulation,  which  was  dictated  by  available  computer 
memory  restrictions.  The  post-processed  bubble  distribution  at  various  cross-section  are 
presented  as  histograms  in  Fig.3.  The  total  bubble  decay  computed  on  the  basis  of  the 
simulations  (Fig.4)  is  consistent  with  the  experimental  data  [8]. 

In  conclusion  we  would  like  to  mention  that  along  with  the  relative  simplicity  of 
implementation,  and  parallel  efficiency  of  the  approach,  this  embarrassingly  parallel  strategy  is 
especially  suitable  for  the  newly  emerging  grid  computing  infrastructure,  which  is  still  not  well 
adapted  for  the  tightly  coupled  problem,  as  those  relying  on  domain  decomposition  methods. 
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Figure  1 :  Particle  ensemble 


Figure  2:  Particle  distribution  in  thb  inlet  plane 
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